這是我研究生課程“現代訊號處理”中的作業報告，上傳到blog中。

---

經典功率譜估計

可以採用直接法，也稱週期圖法，利用公式計算功率譜密度。或者根據自相關函式和譜密度之間的傅立葉變換關係
來計算，稱為間接法或自相關函式法。

還可以先作加窗平滑處，對序列x(n)或估計的自相關函式進行加窗（如漢寧窗、漢明窗）截斷，前者稱作資料窗，後者稱作滯後窗。

---

MATLAB程式設計實現

對訊號x(n)=sin⁡(ωt)+n(t)和x(n)=exp⁡(jωt)+n(t)使用上述4種方法進行功率譜估計。疊加高斯白噪聲n(t)後的訊雜比為-10dB、-3dB、3dB、10dB四種情況。程式設計不使用fft、xcorr等MATLAB提供的函式。

關鍵程式如下，四個函式分別使用不同方法進行功率譜估計，一個函式用於繪圖。呼叫自己編寫的函式實現功能（MATLAB R2017a版本下測試）。

N = 512; fs = 1; M=256;
t = (0:N-1)/fs;
x1n = sin(2*pi*0.1*t);           %訊號1，實訊號
x2n = exp(1j*2*pi*0.1*t);      %訊號2，覆訊號

P_estimation(x1n, 10, N, M, fs);    %估計不同SNR、不同訊號的功率譜
P_estimation(x1n, 3, N, M, fs);      
P_estimation(x1n, -3, N, M, fs);    
P_estimation(x1n, -10, N, M, fs);    
P_estimation(x2n, 10, N, M, fs);    
P_estimation(x2n, 3, N, M, fs);    
P_estimation(x2n, -3, N, M, fs);   
P_estimation(x2n, -10, N, M, fs);    

%   對疊加SNR的訊號xn，使用方法1-4進行功率譜估計
function P_estimation(xn, snr, N, M, fs)
    sn = awgn(xn, snr, 'measured');
    f1=(0:N-1)*fs/N;                        %功率譜座標軸
    Sx1 = direct_method(sn, N);     % 方法1：直接法  
    Sx2 = indir_method(sn,N,M);    % 方法2：間接法  
    Sx3 = add_datawin(sn, N);        % 方法3：資料加窗 
    Sx4 = add_rxwin(sn,N,M);         % 方法4：自相關函式加窗
    figure;
    subplot(221); plot(f1,10*log10(Sx1));
    title('直接法 SNR='+string(snr)); xlabel('f/Hz'); ylabel('Sx(f)(dB/Hz)');
    subplot(222); plot(f1,10*log10(Sx2));
    title('間接法 SNR='+string(snr)); xlabel('f/Hz'); ylabel('Sx(f)(dB/Hz)');
    subplot(223); plot(f1,10*log10(Sx3));
    title('資料加窗 SNR='+string(snr)); xlabel('f/Hz'); ylabel('Sx(f)(dB/Hz)');
    subplot(224); plot(f1,10*log10(Sx4));
    title('自相關函式加窗 SNR='+string(snr)); xlabel('f/Hz'); ylabel('Sx(f)(dB/Hz)');    
end

%    方法1：直接法得到頻率譜估計
function Sx = direct_method(xn, N)   %估計N點序列xn的功率譜Sx
    Sx = zeros(1,N);
    for index = 1 : N
        for n = 0 : N-1       %求Fourier變換
            Sx(index) = Sx(index) + xn(n+1)*exp(-1j*n*index*2*pi/N);  
        end
        Sx(index) = abs(Sx(index))*abs(Sx(index))/N;  %求功率譜
    end
end

%    方法2：間接法得到頻率譜估計
function Sx = indir_method(xn,N,M)   %估計N點序列xn的功率譜Sx    
    Rx = zeros(1,2*M-1);   %M為估計的自相關函式的點數，1<<M<N
    xn = [xn, zeros(1,M)];
    for i = 1 : M                 %估計自相關函式
        for n = 1 : N
            Rx(i+M-1) = Rx(i+M-1) + xn(n+i)*conj(xn(n));
        end
        Rx(i+M-1) = Rx(i+M-1) / N;
    end
    for i = 1 : M-1             %根據共軛關係得到另一半
        Rx(i) = conj(Rx(2*M-i));
    end
    Sx = zeros(1,N);
    for index = 1 : N
        for k = 1 : 2*M-1        %求Fourier變換
            Sx(index) = Sx(index) + Rx(k)*exp(-1j*(k-M)*index*2*pi/N);  
        end
        Sx(index) = abs(Sx(index));
    end
end

%   方法3：資料加窗-修正週期圖
	     在方法1基礎上增加了加窗語句，其餘重複故省略
%  方法4：自相關函式加窗-週期圖平滑
在方法2基礎上增加了加窗語句，其餘重複故省略

---

結果分析

下面是不同訊雜比下對x(n)=sin⁡(ωt)+n(t)訊號的功率譜估計結果。




程式中取樣頻率為1，訊號頻率為0.1。實訊號的功率譜是對稱的，在圖中可以看到兩個尖峰在0.1和與之對稱的0.9處。對於直接法而言，資料加窗後的“尖峰”更窄；對於間接法而言，自相關函式加窗後的功率譜圖更平滑。這也是其稱作“修正週期圖”和“週期圖平滑”的直觀體驗。

下面是不同訊雜比下對x(n)=exp⁡(jωt)+n(t)訊號的功率譜估計結果。




可以得到和上一組測試相同的結論。只不過覆訊號的功率譜不具有對稱性。總體而言，訊雜比越高，訊號的功率峰也就越突出。

---

不同窗函式比較

下面給出直接法和間接法中，不加窗、漢明窗、漢寧窗、布萊克曼窗的差別。為了更直觀地觀察使用不同窗函式之間的區別，沒有新增噪聲。直接法中漢明窗將功率峰變得更窄；漢寧窗和布萊克曼窗效果類似，增大了功率峰值和其餘功率值之間的倍數。

間接法中加三個窗函式都起到了週期圖平滑的作用，漢寧窗和布萊克曼窗的平滑效果看起來比漢明窗要好一些。但三個窗函式都沒有改變功率峰值和其餘功率值之間的倍數關係。