頻域特徵值提取的MATLAB代碼實現（頻譜、功率譜、倒頻譜）

來自專欄 與信號處理有關的那些東東

本篇文章是信號的各種頻域分析方法的理解（頻譜、能量譜、功率譜、倒頻譜、小波分析的後續，對文章中提到的頻域分析方法進行代碼實現。

時域特徵值的代碼實現可以參考：時域特徵值提取的MATLAB代碼實現（均方根、峰值因子、脈衝因子、裕度因子、峭度因子、波形因子和偏度等）

文章如要轉載請私信與我聯繫，並註明來源知乎專欄與信號處理有關的那些東東作者Mr.括弧。

一、頻譜

頻譜用到的函數主要是fft和fftshift。

需要注意的主要有三點：

1.直接做fft的結果，信號的前半部分對應頻率[0,fs/2]，後半部分對應[-fs/2,0]。參見頻譜結果圖的第2張。為了將零頻點移到頻譜中間，需要使用fftshift函數，結果參見頻譜結果圖的第3張。

2.通常我們關心的都是正頻率區間的結果，有兩種截取方法，一種是在fftshift的結果中截後半段，一種是在fft的結果中截前半段，其結果是一樣的。後一種方法更簡潔。具體參見頻譜結果圖的第4、5張。

3.根據奈奎斯特定理，信號的採樣頻率（1/t_s）必須大於信號頻率最大值的兩倍。

t_s = 0.01; %採樣周期t_start = 0.5; %起始時間t_end = 5; %結束時間t = t_start : t_s : t_end;y = 1.5*sin(2*pi*5*t)+3*sin(2*pi*20*t)+randn(1,length(t)); %生成信號%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%頻譜%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%y_f = fft(y); %傅里葉變換subplot(5,1,1);plot(t,y);title(original signal); %繪製原始信號圖Druation = t_end -t_start; %計算採樣時間Sampling_points = Druation/t_s +1; %採樣點數，fft後的點數就是這個數f_s = 1/t_s; %採樣頻率f_x = 0:f_s/(Sampling_points -1):f_s; %注意這裡和橫坐標頻率對應上了，頻率解析度就是f_s/(Sampling_points -1)t2 = f_x-f_s/2;shift_f = abs(fftshift(y_f));subplot(5,1,2);plot(f_x,abs(y_f));title(fft transform);subplot(5,1,3);plot(f_x-f_s/2,shift_f);title(shift fft transform); %將0頻率分量移到坐標中心subplot(5,1,4);plot(t2(length(t2)/2:length(t2)),shift_f(length(shift_f)/2:length(shift_f)));title(shift fft transform); %保留正頻率部分subplot(5,1,5);plot(f_x(1:length(f_x)/2),abs(y_f(1:length(f_x)/2)));title(fft cut); %直接截取fft結果的前半部分

二、功率譜

功率譜有兩種求法：1.(傅立葉變換的平方)/(區間長度)；2.自相關函數的傅里葉變換。這兩種方法分別叫做直接法和相關函數法。（參見信號的各種頻域分析方法的理解（頻譜、能量譜、功率譜、倒頻譜、小波分析））

下述代碼中，直接法就用了(傅立葉變換的平方)/(區間長度)的方法求解的，其結果和使用MATLAB的函數periodogram（周期圖法）結果相同；雖然理論上直接法和相關函數法相同，不過模擬結果中相關函數法對雜訊的抑制效果更好，圖線更平滑。

Fs = 1000;nfft = 1000; %fft採樣點數%產生序列n = 0:1/Fs:1;xn = cos(2*pi*100*n) + 3*cos(2*pi*200*n)+(randn(size(n)));subplot(5,1,1);plot(xn);title(加噪信號);xlim([0 1000]);grid on%FFTY = fft(xn,nfft);Y = abs(Y);subplot(5,1,2);plot((10*log10(Y(1:nfft/2))));title(FFT);xlim([0 500]);grid on%FFT直接平方Y2 = Y.^2/(nfft);subplot(5,1,3);plot(10*log10(Y2(1:nfft/2)));title(直接法);xlim([0 500]);grid on%周期圖法window = boxcar(length(xn)); %矩形窗[psd1,f] = periodogram(xn,window,nfft,Fs);psd1 = psd1 / max(psd1);subplot(5,1,4);plot(f,10*log10(psd1));title(周期圖法);ylim([-60 10]);grid on%自相關結果cxn = xcorr(xn,unbiased); %計算自相關函數%自相關法CXk = fft(cxn,nfft);psd2 = abs(CXk);index = 0:round(nfft/2-1);k = index*Fs/nfft;psd2 = psd2/max(psd2);psd2 = 10*log10(psd2(index+1));subplot(5,1,5);plot(k,psd2);title(間接法);grid on

下圖中的縱坐標都進行了取對數的處理（10log），取對數的目的是使那些振幅較低的成分相對高振幅成分得以拉高，以便觀察掩蓋在低幅雜訊中的信號特徵。

三、倒頻譜

倒頻譜的求解函數為rceps（實倒頻譜），在MATLAB的幫助文檔中，rceps的計算公式為real(ifft(log(abs(fft(y)))))，即信號→頻譜→對數→傅里葉逆變換，而倒頻譜的定義表述中卻是信號→功率譜→對數→傅里葉逆變換。即功率譜被換成了頻譜。私以為是因為功率譜為頻譜值的平方，在取對數後平方會變成係數2，對後續計算影響不大，因而可以近似認為結果相同。

在模擬中要看出倒頻譜的作用，需要手動生成一組調製信號。下列程序將高頻（主頻為50/100/200Hz）和低頻（主頻為5/10/20Hz）信號進行調製，分別畫出低頻、高頻和調製信號的時域圖和頻譜圖。在圖ftt_y中可以看到邊頻帶的形成。（邊頻帶相關概念參見信號的各種頻域分析方法的理解（頻譜、能量譜、功率譜、倒頻譜、小波分析）

代碼如下：

sf = 1000;nfft = 1000;x = 0:1/sf:5;y1=10*cos(2*pi*5*x)+7*cos(2*pi*10*x)+5*cos(2*pi*20*x)+0.5*randn(size(x));y2=20*cos(2*pi*50*x)+15*cos(2*pi*100*x)+25*cos(2*pi*200*x)+0.5*randn(size(x));for i = 1:length(x) y(i) = y1(i)*y2(i);endsubplot(3,3,1)plot(y1);xlim([0 5000]);title(y1);subplot(3,3,2)plot(y2);xlim([0 5000]);title(y2);subplot(3,3,3)plot(y);xlim([0 5000]);title(y=y1*y2);t = 0:1/sfnfft-1)/sf;nn = 1:nfft;subplot(3,3,4)ft = fft(y1,nfft);Y = abs(ft);plot(0:nfft/2-1,((Y(1:nfft/2))));title(fft_y_1);ylabel(幅值);xlim([0 300]);grid on;subplot(3,3,5)ft = fft(y2,nfft);Y = abs(ft);plot(0:nfft/2-1,((Y(1:nfft/2))));title(fft_y_2);ylabel(幅值);xlim([0 300]);grid on;subplot(3,3,6)ft = fft(y,nfft);Y = abs(ft);plot(0:nfft/2-1,((Y(1:nfft/2))));title(fft_y);ylabel(幅值);xlim([0 300]);grid on;subplot(3,3,7)z = rceps(y);plot(t(nn),abs(z(nn)));title(z=rceps(y));ylim([0 0.3]);xlabel(時間(s));ylabel(幅值);grid on;subplot(3,3,8)yy = real(ifft(log(abs(fft(y))))); %信號→傅里葉→對數→傅里葉逆變換plot(t(nn),abs(yy(nn)));title(real(ifft(log(abs(fft(y))))));ylim([0 0.3]);xlabel(時間(s));ylabel(幅值);grid on;

上圖紅圈處可以看到三個峰值，分別是0.05s、0.1s和0.2s。其對應的頻率分別為20Hz、10Hz和5Hz，正是調製信號的低頻分量。該低頻分量在圖fft_y中以邊頻帶的形式出現，無法看出其對應頻率值，而在倒頻譜中能輕易展現。這正是倒頻譜的意義所在。

參考：

對功率譜的一點理解 - cynchanpin - 博客園

信號處理方法 - 聲振論壇 - 振動,動力學,聲學,信號處理,故障診斷 - Powered by Discuz!

推薦閱讀：