clear
numberOfPoints = 1024;                           % 采樣點(diǎn)數(shù)
fs = 44100;                                      % 采樣頻率  
fsOfSignal = 500;                                % 原信號(hào)的估計(jì)頻率
ch = 1;                                          % 聲音采樣聲道數(shù)

maxNumber = (1/fs) * (numberOfPoints - 1);
maxNumber2 = (1/fs) * (numberOfPoints/2 - 1);
x  = 0 : 1/fs : maxNumber;
x2 = 0 : 1/fs : maxNumber2;                      % 產(chǎn)生時(shí)域橫軸


% RANDOM SIGNAL

% wave = wavrecord(numberOfPoints * 2,fs,ch,'double');
% for i = 1:1024
%     waveSignal(i) = wave(i+1024);
% end                                            % 聲音信號(hào)采樣

waveSignal = sin(fsOfSignal * 2 * pi * x);       % 用來(lái)模擬的周期正弦信號(hào)

ranSignal  = awgn(waveSignal,-7);                % 為周期信號(hào)疊加高斯白噪聲 信噪比 1:2


figure(1)                                        % 顯示原周期信號(hào)和混合高斯白噪聲后的信號(hào)

subplot(2,1,1);
plot(x,waveSignal,'b-');
grid on;
title('原周期信號(hào)')

subplot(2,1,2);
plot(x,ranSignal,'b-');
grid on;
title('混合信號(hào)')

haveSignal = 1 - textWgn(ranSignal,numberOfPoints/2,numberOfPoints)    % 檢測(cè)混合信號(hào)中是否存在有用信號(hào) 1表示有 0表示沒有

equRS   = equality(ranSignal,numberOfPoints)                           
mVarRS  = meanVar(ranSignal,numberOfPoints)                           
varRS   = variance(equRS,mVarRS)                                      % 計(jì)算混合信號(hào)的均值、均方值、方差

figure(2);                                                             % 顯示混合信號(hào)的均值、方差、頻譜、自相關(guān)函數(shù)、功率譜

subplot(4,1,1);
plot(x,equRS,'b-');
grid on;
title('混合信號(hào)均值')

subplot(4,1,2);
plot(x,varRS,'b-');
grid on;
title('混合信號(hào)方差')

frequencyRS = fft(ranSignal);                                         % 計(jì)算混合信號(hào)的1024點(diǎn)FFT變換

dfRS = (1/(maxNumber/(numberOfPoints-1)))/numberOfPoints;
nRS  = 0:numberOfPoints-1;
fRS  = nRS*dfRS;                                                      % 進(jìn)行橫坐標(biāo)的頻域變換

subplot(4,1,3);
plot(fRS,abs(frequencyRS(nRS+1))*2/numberOfPoints,'r-');              % 進(jìn)行縱坐標(biāo)的修正，畫出頻譜
grid on;
title('混合信號(hào)頻譜');
axis([0,2.25e4,0,1])

selfRelationRS = relation(ranSignal,ranSignal,equRS,equRS,varRS,varRS,numberOfPoints/2,numberOfPoints);  % 計(jì)算自相關(guān)函數(shù)

powerRS = fft(selfRelationRS);                                        % 功率譜密度是自相關(guān)函數(shù)的傅立葉變換
pdfRS   = (1/(maxNumber/(numberOfPoints/2-1)))/(numberOfPoints/2);
pnRS    = 0:numberOfPoints/2-1;
pfRS    = pnRS*pdfRS;

subplot(4,1,4);
plot(pfRS,abs(powerRS(pnRS+1))*2/(numberOfPoints/2),'k-');            % 計(jì)算功率譜密度
grid on;
title('功率譜密度');
axis([0,1.25e4,0,0.2])

figure(3);                                                            % 顯示經(jīng)過濾波后的信號(hào)，與原信號(hào)進(jìn)行對(duì)比

equWS   = equality(waveSignal,numberOfPoints)
mVarWS  = meanVar(waveSignal,numberOfPoints)
varWS   = variance(equWS,mVarWS)
selfRelationWS = relation(waveSignal,waveSignal,equWS,equWS,varWS,varWS,numberOfPoints/2,numberOfPoints);   % 計(jì)算原信號(hào)各項(xiàng)數(shù)字特征

[b,a] = butter(9,fsOfSignal / (fs/2),'low');                           % Butterworth 低通濾波器 3dB截止頻率等于原信號(hào)估計(jì)頻率
freqz(b,a,128,1000);                                                  % 顯示濾波器的幅頻特性和相頻特性

signalRE = filter(b,a,ranSignal);                                     % 進(jìn)行濾波，恢復(fù)原信號(hào)

equRE   = equality(signalRE,numberOfPoints)                        
mVarRE  = meanVar(signalRE,numberOfPoints)
varRE   = variance(equRE,mVarRE)
selfRelationRE = relation(signalRE,signalRE,equRE,equRE,varRE,varRE,numberOfPoints/2,numberOfPoints);     % 計(jì)算恢復(fù)信號(hào)的各項(xiàng)數(shù)字特性

figure(4)

subplot(2,1,1)
plot(x,waveSignal,'b-');
grid on;
title('原信號(hào)')

subplot(2,1,2);
plot(x,signalRE,'b-');
grid on;
title('經(jīng)過濾波后的信號(hào)')

figure(5);                                                         % 顯示原信號(hào)各項(xiàng)數(shù)字特征

subplot(5,1,1);
plot(x,equWS,'b-');
grid on;
title('原信號(hào)均值')

subplot(5,1,2);
plot(x,varWS,'b-');
grid on;
title('原信號(hào)方差')

frequencyWS = fft(waveSignal);

dfWS = (1/(maxNumber/(numberOfPoints-1)))/numberOfPoints;
nWS  = 0:numberOfPoints-1;
fWS  = nWS*dfWS;
selfRelationWS = relation(waveSignal,waveSignal,equWS,equWS,varWS,varWS,numberOfPoints/2,numberOfPoints);

subplot(5,1,3);
plot(x2,selfRelationWS)
grid on;
title('原信號(hào)自相關(guān)函數(shù)')

subplot(5,1,4);
plot(fWS,abs(frequencyWS(nWS+1))*2/numberOfPoints,'r-');
grid on;
title('原信號(hào)頻譜');
axis([0,2.25e4,0,1]);

powerWS = fft(selfRelationWS);
pdfWS   = (1/(maxNumber/(numberOfPoints/2-1)))/(numberOfPoints/2);
pnWS    = 0:numberOfPoints/2-1;
pfWS    = pnWS*pdfWS;

subplot(5,1,5);
plot(pfWS,abs(powerWS(pnWS+1))*2/(numberOfPoints/2),'k-');
grid on;
title('原信號(hào)功率譜密度');
axis([0,1.25e4,0,1]);

figure(6);                                                       % 顯示恢復(fù)信號(hào)的各項(xiàng)數(shù)字特征

subplot(5,1,1);
plot(x,equRE,'b-');
grid on;
title('經(jīng)過濾波后信號(hào)的均值')

subplot(5,1,2);
plot(x,varRE,'b-');
grid on;
title('經(jīng)過濾波后信號(hào)的方差')

frequencyRE = fft(signalRE);

dfRE = (1/(maxNumber/(numberOfPoints-1)))/numberOfPoints;
nRE  = 0:numberOfPoints-1;
fRE  = nRE*dfRE;
selfRelationRE = relation(signalRE,signalRE,equRE,equRE,varRE,varRE,numberOfPoints/2,numberOfPoints);

subplot(5,1,3);
plot(x2,selfRelationRE);
grid on;
title('經(jīng)過濾波后信號(hào)的自相關(guān)函數(shù)')

subplot(5,1,4);
plot(fRE,abs(frequencyRE(nRE+1))*2/numberOfPoints,'r-');
grid on;
title('經(jīng)過濾波后信號(hào)的頻譜');
axis([0,2.25e4,0,0.5]);


powerRE = fft(selfRelationRE);
pdfRE   = (1/(maxNumber/(numberOfPoints/2-1)))/(numberOfPoints/2);
pnRE    = 0:numberOfPoints/2-1;
pfRE    = pnRE*pdfRE;

subplot(5,1,5);
plot(pfRE,abs(powerRE(pnRE+1))*2/(numberOfPoints/2),'k-');
grid on;
axis([0,1.25e4,0,2]);

title('經(jīng)過濾波后信號(hào)的功率譜密度');
 
periodWS = 1/fsOfSignal                                            % 計(jì)算原信號(hào)周期
[maxFsRE,i] = max(abs(frequencyRE(nRE+1))*2/numberOfPoints);
periodRE = 1/((i/numberOfPoints)*fs)                               % 計(jì)算恢復(fù)信號(hào)周期