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('自相關函數信號圖形');