%信號的傳送，提取及分析
f=1000;%正弦信號頻率
F=44000;%采樣頻率
A=1;
n=1:F;
s=wavrecord(44000,44000,1);
s=s';
%s=A*cos(2*pi*f*n/F);%信號產生

%S=fftshift(fft(s));

%plot(2/F*n-1,abs(S),'*')

figure(1)
plot(n,s);%繪制信號波形
axis([1,3*F/f,-2*A,2*A]);
title('signal'),xlabel('t'),ylabel('s(t)');

sn=s+normrnd(0,A,1,F);%信號＋噪聲

figure(2)
t=1:fix(3*F/f);
psn=sn(1:3*F/f);
plot(t,psn);%繪制信號＋噪聲波形
title('signal+noise'),xlabel('t'),ylabel('s(t)+n(t)');

Y=fftshift(fft(sn));%信號＋噪聲的頻譜

figure(3)
plot(2*(n-1)/F-1,abs(Y));%繪制信號＋噪聲的頻譜
title('signal+noise fourier transform'),xlabel('w/pi'),ylabel('S(w)+N(w)');

p=5000;%低通濾波器帶寬
f1=[0,2*p/F,2*p/F,1];
m=[1 1 0 0];
b=fir2(40,f1,m);%低通濾波器設計

figure(4)%幅頻響應
[h,w]=freqz(b,1,128);
plot(w/pi,abs(h));
title('lowpass filter'),xlabel('w/pi'),ylabel('H(w)');

sn1=conv(sn,b);%信號通過低通濾波器
Y1=fftshift(fft(sn1));%信號＋噪聲通過低通濾波器的頻譜

figure(5)
t=1:max(size(Y1));
plot(2/max(size(Y1))*t-1,abs(Y1));%繪制信號＋噪聲通過低通濾波器的頻譜
title('signal+noise passing the lowpass filter'),xlabel('w/pi'),ylabel('S(w)+N(w)');

for i=1:fix(F/f)                   %平均法求得信號
    sum=0;
    for j=1:f
    sum=sum+sn1(i+fix(F/f)*(j-1));
    end
    signal(i)=sum/f;
end

sign=[signal signal signal];
figure(6)
plot(sign)
S=fftshift(fft(signal));
figure(7)
t=1:max(size(S));
plot(2*(t-1)/(F/f)-1,abs(S));
title('signal recovery'),ylabel('s1(t)');

disp('信號均值');
Ex1=mean(signal) 
disp('信號方差值');
D=var(signal)
%信號和噪聲的自相關函數
disp('信號和噪聲自相關函數');
R=xcorr(signal,signal);


figure(8)%自相關函數圖像
t=1:max(size(R));
plot(t-44,R);
%axis([1000,2000,-0.6,0.6])
title('自相關函數示意圖像');