monew11.m

给定一个音频范围内的模拟强噪声随机信号（强噪声+周期信号）

clear all;close all;
N=4096;    %采样点个数         %获得信号，并对其频谱自相关函数进行分析
Fs=44100;   %采样频率
f=1000;   %信号频率
M = fix(Fs/f);   %一个周期内采样的点数
n=1:1:N;   %采样向量
n=n';   %将行向量转置成列向量
t=n/Fs;         %采样时间

y =20*wavrecord(N,Fs,1,'double');%从声卡采样初始信号          
%y =10*sin(2*f*pi*t);
figure(1);
subplot(311);
plot(t,y);axis([0,0.01,-11,11]);         
xlabel('原始信号   y(t)  /0.1s');
grid
yk=fftshift(fft(y,N));                                  
k=-(N-1)/2:(N-1)/2;   %求对称于零频率的FFT位置向量
f1=k*Fs/N;%求幅度频率中的频率向量
subplot(312);
plot(f1,abs(yk));axis([0,5000,0,3*10^4]);
xlabel('原始信号频谱   y(f)   /Hz');
[C0,LAGS0] = xcorr(y,'unbiased');%求自相关函数   无偏估计的互相关函数

subplot(313);
plot(LAGS0/Fs,C0);axis([0,0.01,-100,100]);
xlabel('自相关函数    t /s');ylabel('Ry(t)');
grid

fprintf('原始信号的均值为：');
y_average = mean(y)
fprintf('原始信号的方差为：');
y_var = var(y)
fprintf('原始信号的平均功率为：');
y_Power_average = sum(abs(y).^2)/length(y) 




z =awgn(y,-20);                     %加入噪声
                  %合成信号，分析其特性
                  
s=10;   %取前10个周期的长度来恢复原信号    %通过求多个周期内的和的均值恢复信号
D=s*M;        %返回原信号时原信号的点数
K = N/M;      %取样的总周期
zz=z(1:D,1);  %取前10周期内所有点
for a=1:(K-1-s)
    if a*M+D<N
    za=z((a*M+1):(a*M+D),1);
    zz=zz+za;%累加各个周期的对应点
    end
end
yy=zz/(K-s);%求各采样点的均值,从而消除高斯噪声

n=1:1:D;
t=n/Fs;
figure(2);
subplot(311);
plot(t,yy);axis([0,0.011,-10,10]);
xlabel('求平均后的信号波形  yy(t)   /s');
yyk=fftshift(fft(yy));
k=-(D-1)/2:(D-1)/2;   %求对称于零频率的FFT位置向量
f3=k*Fs/D;        %求幅度频率中的频率向量
subplot(312);
plot(f3,abs(yyk));axis([0,5000,0,3000]);
xlabel('求平均后信号的频谱特性');
[C3,LAGS3] = xcorr(yy,'unbiased');%求自相关函数
subplot(313);
plot(LAGS3/Fs,C3);axis([0,0.01,-100,100]);
xlabel('求平均后信号的自相关函数    t /s');ylabel('Ryy(t)');                  
 grid                 
                  
                  
                  
n=1:1:N;%采样向量
n=n';%将行向量转置成列向量
t =n/Fs; %采样时间
figure(3);
subplot(411);
plot(t,z);
xlabel('合成信号   z(t)  /s');
subplot(412);
plot(t,z);axis([0 0.01 -10,10]);
xlabel('合成信号   z(t)  /0.1s');
zk=fftshift(fft(z,N));
k=-(N-1)/2:(N-1)/2;   %求对称于零频率的FFT位置向量
f2=k*Fs/N;        %求幅度频率中的频率向量
subplot(413);
plot(f2,abs(zk));axis([0,5000,0,3*10^4]);
xlabel('合成信号频谱   z(w)   /Hz');
[C2,LAGS2] = xcorr(z,'unbiased');%求自相关函数
subplot(414);
plot(LAGS2/Fs,C2);
xlabel('合成信号的自相关函数    t /s');ylabel('Rz(t)');
axis([0,0.01,-100,100]);
grid

fprintf('合成信号的均值为：');
z_average = mean(z)
fprintf('合成信号的方差为：');
z_var = var(z)
fprintf('合成信号的平均功率为：');
z_Power_average = sum(abs(z).^2)/length(z) 


Rp=0.1;Rs=80;
f=[1100,1500];m=[1,0];
dev=[(10^(Rp/20)-1)/(10^(Rp/20)+1),10^(-Rs/20)];
[N,fo,mo,w]=remezord(f,m,dev,Fs);
b=remez(N,fo,mo,w);
[hf,f]=freqz(b,1,4096,Fs);
figure(4);
subplot(511);
plot(f,abs(hf));axis([0,5000,0,1.2]);
xlabel('数字滤波器频谱特性');
ht=ifft(hf);
yyy=conv(z,ht);
n=1:1:length(yyy);   %采样向量
n=n';      %将行向量转置成列向量
t =n/Fs;         %采样时间
subplot(512);
plot(t,2*yyy);axis([0,0.1,-10,10]);
xlabel('通过数字滤波器后的信号');


subplot(513);
plot(f1,abs(yk));axis([0,5000,0,3*10^4]);

xlabel('通过数字滤波器后的信号频谱   yyy(f)   /Hz');

[C4,LAGS4]=xcorr(yyy,'unbiased'); 

f=1000;
subplot(514);
g=4*C4;
plot(LAGS4/Fs,g);axis([0,0.01,-100,100]);
xlabel('通过滤波器后的相关函数');
grid

v=max(g)
A=(v*2)^0.5

t=n/Fs;
f=1000;
k =A*sin(2*f*pi*t);
subplot(515);
plot(t,k);axis([0,0.01,-11,11]);
xlabel('由相关函数恢复出的信号');
grid