m=80000;
n=44100;
x=0:1/n:79999/n;
Signal=wavrecord(m,n,'double');
Signal=Signal';

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

F1=fftshift(fft(Signal));%頻譜函數(shù)數(shù)值
figure(2);
g = floor((-length(F1)/2+0.5):(length(F1)/2-0.5));
d=g*n/length(F1);
plot(d,abs(F1));
axis([0,10000,0,100]);
title('信號頻譜');

GNoise=wgn(1,m,-20);
figure(3)
plot(x,GNoise);
axis([0,0.015,-2,2]);
title('高斯噪聲');

F2=fftshift(fft(GNoise));%頻譜函數(shù)數(shù)值
figure(4);
g = floor((-length(F2)/2+0.5):(length(F2)/2-0.5));
d=g*44100/length(F2);
plot(d,abs(F2));
axis([0,10000,0,100]);
title('噪聲頻譜');


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

F3=fftshift(fft(V));%頻譜函數(shù)數(shù)值
figure(6);
g = floor((-length(F3)/2+0.5):(length(F3)/2-0.5));
d=g*n/length(F3);
plot(d,abs(F3));
axis([0,10000,0,100]);
title('信號+噪聲 頻譜');

disp('自相關(guān)函數(shù)');
lv=xcorr(V,V,'unbiased');
figure(7);
plot(lv);axis([1000,1500,-0.3,0.3])
title('自相關(guān)函數(shù)示意圖像');

%figure(20)
%lv3=xcorr(Signal,Signal,'unbiased')
%plot(lv3);axis([1000,2000,-0.3,0.3])
%title('測試信號');

figure(18);
lv1=10*diff(lv);
plot(lv1);
axis([1000,1500,-0.2,0.2]);
title('恢復(fù)信號');


F4=fftshift(fft(lv));%頻譜函數(shù)數(shù)值
figure(8);
g = floor((-length(F4)/2+0.5):(length(F4)/2-0.5));
d=g*n/length(F4);
plot(d,abs(F4));
axis([0,10000,0,100]);
title('自相關(guān)函數(shù)恢復(fù)信號頻譜');

%各函數(shù)功率譜
window=hamming(500);
noverlap=20;
range='half';
fs=n;

[Py,f]=pwelch(Signal,window,noverlap,m,fs,range);
[Py1,f]=pwelch(GNoise,window,noverlap,m,fs,range);
[Py2,f]=pwelch(V,window,noverlap,m,fs,range);
[Py3,f]=pwelch(lv,window,noverlap,m,fs,range);

DBPy=10*log10(Py);
DBPy1=10*log10(Py1);
DBPy2=10*log10(Py2);
DBPy3=10*log10(Py3);

figure(9)
plot(f,DBPy);
xlabel('F(Hz)');
ylabel('DB');
title('原信號功率譜');

figure(10)
plot(f,DBPy1);
xlabel('F(Hz)');
ylabel('DB');
title('噪聲信號功率譜');

figure(11)
plot(f,DBPy2);
xlabel('F(Hz)');
ylabel('DB');
title('信號+噪聲功率譜');

figure(12);
plot(f,DBPy3);
xlabel('F(Hz)');
ylabel('DB');
title('恢復(fù)信號功率譜');


%信號分析

%信號均值
disp('信號均值');
Ex1=mean(lv) 
%信號方差值
disp('信號方差值');
D=var(lv)
%信號標(biāo)準(zhǔn)差
disp('信號標(biāo)準(zhǔn)差');
S=std(lv)
%信號均方值
disp('信號均方值');
Ex2=Ex1^2+D



%恢復(fù)信號分析

%信號均值
disp('信號均值');
Ex1=mean(lv) 
%信號方差值
disp('信號方差值');
D=var(lv)
%信號標(biāo)準(zhǔn)差
disp('信號標(biāo)準(zhǔn)差');
S=std(lv)
%信號均方值
disp('信號均方值');
Ex2=Ex1^2+D