niub.m

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

x=0:1/44100:2047/44100;
Signal=wavrecord(2048,44100,1,'double');
Signal=Signal';

figure(1)
plot(x,Signal);
axis([0,0.015,-1,1]);
title('源信号');


Fx=fftshift(fft(Signal));%频谱函数数值
figure(2);
g = floor((-length(Fx)/2+0.5):(length(Fx)/2-0.5));
% D=2*pi/length(Fx);
% g=g*D*fs/2*pi;
d=g*44100/length(Fx);
plot(d,abs(Fx));
axis([-3000,3000,0,200]);
xlabel('输入信号频谱函数');



wp=1000*2*pi;ws=2000*2*pi;Fs=44100;Rp=0.3;Rs=50;
[N,Wn]=cheb1ord(wp,ws,Rp,Rs,'s');
[z,p,k]=cheb1ap(N,Rp);
[A,B,C,D]=zp2ss(z,p,k);
[AT,BT,CT,DT]=lp2lp(A,B,C,D,Wn);
[b,a]=ss2tf(AT,BT,CT,DT);
[bz,az]=impinvar(b,a,Fs);
[H,W]=freqz(bz,az);
figure(3);
plot(W*Fs/(2*pi),abs(H));grid;
title('滤波器的幅频特性');
xlabel('频率/Hz');
ylabel('幅度');
axis([0 2500 0 1.2]);
lv=filter(bz,az,Signal);
figure(4);
plot(lv);
axis([200,800,-1,1]);
title('通过低通滤波器后的信号');

GNoise=wgn(1,2048,-10);
figure(5)
plot(x,GNoise);
axis([0,0.015,-1,1]);
title('高斯噪声');

Fx=fftshift(fft(GNoise));%频谱函数数值
figure(6);
g = floor((-length(Fx)/2+0.5):(length(Fx)/2-0.5));
% D=2*pi/length(Fx);
% g=g*D*fs/2*pi;
d=g*44100/length(Fx);
plot(d,abs(Fx));
axis([-3000,3000,0,100]);
xlabel('高斯噪声频谱函数');

V=Signal+GNoise;
figure(7)
plot(x,V);
axis([0,0.015,-1,1]);
title('加噪后的信号');

sin=filter(bz,az,V);
figure(8)
plot(x,sin,'r');
hold on;
plot(x,V);
axis([0,0.015,-2,2]);
title('滤波后信号与带噪信号的比较,红线为滤波后信号');
hold off;
Fx=fftshift(fft(sin));%频谱函数数值
figure(9);
g = floor((-length(Fx)/2+0.5):(length(Fx)/2-0.5));
% D=2*pi/length(Fx);
% g=g*D*fs/2*pi;
d=g*44100/length(Fx);
plot(d,abs(Fx));
axis([-5000,5000,0,200]);
xlabel('带噪信号滤波后频谱函数');


Fx=fftshift(fft(V));%频谱函数数值
figure(10);
g = floor((-length(Fx)/2+0.5):(length(Fx)/2-0.5));
% D=2*pi/length(Fx);
% g=g*D*fs/2*pi;
d=g*44100/length(Fx);
plot(d,abs(Fx));
axis([-5000,5000,0,200]);
xlabel('带噪信号频谱函数');

%计算信号的均值、方差、
%信号均值
disp('信号均值');
Ex1=mean(lv) 
%信号方差值
disp('信号方差值');
D=var(lv)
%信号标准差
disp('信号标准差');
S=std(lv)
%信号均方值
disp('信号均方值');
Ex2=Ex1^2+D
%信号和噪声的自相关函数
disp('信号和噪声自相关函数');
R=xcorr(V,V,'unbiased');


%分析处理信号的频谱
disp('频谱函数数值');
Fx=fftshift(fft(R));%频谱函数数值
figure(12);
g = floor((-length(Fx)/2+0.5):(length(Fx)/2-0.5));
% D=2*pi/length(Fx);
% g=g*D*fs/2*pi;
d=g*44100/length(Fx);
plot(d,abs(Fx));
axis([-3000,3000,0,200]);
xlabel('FFt频谱函数');

figure(13);
plot(R);
axis([200,800,-0.1,0.1]);
title('自相关函数信号图形');