求功率譜有好幾種方法,本文例兩種,一種是先求訊號頻譜,然後取平方後平均;第二種是用welch法(可指定各種窗)。
rng default;
fs = 1000; % 取樣率
t = 0:1/fs:1-1/fs;
% x = cos(2*pi*50*t)+ randn(size(t));
x = cos(2*pi*50*t);
nfft = length(x); %指定傅利葉變換的點數(點數越大頻率解析度越高)
figure(1); %畫單邊頻譜圖
sig_dft = abs(fft(x, nfft));
sig_dft(2:end-1) = 2*sig_dft(2:end-1); %因為只顯示單邊譜,所以需要把數值*2,摺起來
sig_amp = mapminmax(sig_dft/nfft); %歸一化
x1 = (0:nfft-1)*fs/nfft;
plot(x1(1:250),sig_amp(1:250));
xlabel('frequency(hz)');
ylabel('amplitude');
title('single sided amplitude spectrum');
figure(2); %常規方法畫單邊功率譜
xdft = fft(x);
xdft = xdft(1:nfft/2+1);
psdx = (1/(fs*nfft)) * abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);
freq = 0:fs/nfft:fs/2;
semilogy(freq,psdx); %semilogy以對數方式繪製縱軸
figure(3); %修正週期圖功率譜密度估計方法
[p,f]=pwelch(x, hanning(nfft), nfft/2, nfft, fs); %指定不同的窗會影響功率譜結果
semilogy(f, p);
xlabel('frequency(hz)');
ylabel('power spectrum density (db/hz)');
title('single sided power spectrum density');
figure(4); %反變換,顯示原影象
h = ifft(sig_dft, nfft);
plot(real(h)); %指定顯示實部
按道理講,figure2 和figure3的影象理應一致。但是在figure 3 這裡我並不知道在welch法中應該指定什麼窗才能得到figure 2一樣的結果…有大佬知道的還望賜教
[1]matlab求功率譜密度(psd)官方示例 (無法跳轉請將**貼上到搜尋欄)
[2]matlab功率譜函式psd與pwelch等價問題
[3]matlab 求功率譜: 功率譜密度的數學定義
[4]頻譜縱軸、單邊譜部分解釋
[6]知乎:功率譜密度函式
能量譜與功率譜
能量譜與功率譜分別是針對能量有限的訊號和功率有限的訊號。在進行訊號的譜分析的時候一定更要看準這是乙個能量訊號還是乙個功率訊號,應用不同的譜進行分析會使問題的解決思路更加明確。對於週期功率訊號我們使用傅利葉級數的方法進行分析。指數形式的傅利葉級數係數的求解用下面的公式。對於非週期的能量訊號,我們主要使...
功率譜密度
乙個訊號x t 是乙個波,這個波必然有能量。而將這個波分解成很多個波以後,每個波都乙個能量,所有這些波的能量相加的和應是原來訊號的能量。如將乙個訊號x t 用傅利葉變換分解成頻域。可以認為x t 是乙個週期是 00到 00的函式,有很多種分解方法,如可以將其分解為所有頻率是一維實數集的形式。而正常的...
頻譜 幅度譜 功率譜和能量譜
在訊號處理的學習中,有一些與譜有關的概念,如頻譜 幅度譜 功率譜和能量譜等,常常讓人很糊塗,搞不清其中的關係。這裡主要從概念上釐清其間的區別。對乙個時域訊號進行傅利葉變換,就可以得到的訊號的頻譜,訊號的頻譜由兩部分構成 幅度譜和相位譜。這個關係倒還是簡單。那麼,什麼是功率譜呢?什麼又是能量譜呢?功率...