在信號(hào)處理領(lǐng)域，頻域分析是至關(guān)重要的，但很多人在編程實(shí)現(xiàn)這一過程中遇到了困難。這恰恰是我們今天需要重點(diǎn)探討的核心議題。

直接FFT的結(jié)果與調(diào)整

在執(zhí)行FFT分析后，結(jié)果呈現(xiàn)出了明顯的規(guī)律。例如，信號(hào)的起始部分對(duì)應(yīng)的是頻率區(qū)間[0,fs/2]，而結(jié)束部分則對(duì)應(yīng)[-fs/2,0]。在實(shí)際操作中，必須將零頻點(diǎn)放置在頻譜的中央，這非得借助fftshift函數(shù)不可。以一個(gè)具體項(xiàng)目為例，在分析一段采集到的音頻信號(hào)時(shí)，未經(jīng)fftshift處理的FFT結(jié)果讓人難以直接識(shí)別頻率分布。但經(jīng)過調(diào)整，頻譜圖便清晰地展現(xiàn)了正負(fù)頻率的分布。這一步驟對(duì)于聲音信號(hào)的頻域分析等類似工作極為關(guān)鍵。

我們得知，處理FFT結(jié)果需要采用截取正頻段的方法。這主要有兩種方式：一是在fftshift后取后半部分，二是直接在fft結(jié)果中取前半部分。這兩種方法得出的結(jié)果并無差異。就好比分析心電圖信號(hào)，不管選用哪種截取方式，都是為了提取正頻段成分，以便進(jìn)行與疾病相關(guān)的頻率特征分析。

t_s = 0.01; %采樣周期
t_start = 0.5; %起始時(shí)間
t_end = 5;     %結(jié)束時(shí)間
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));  %生成信號(hào)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%頻譜%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y_f = fft(y); %傅里葉變換
subplot(5,1,1);
plot(t,y);title('original signal');   %繪制原始信號(hào)圖
Druation = t_end -t_start;  %計(jì)算采樣時(shí)間
Sampling_points = Druation/t_s +1;  %采樣點(diǎn)數(shù)，fft后的點(diǎn)數(shù)就是這個(gè)數(shù)
f_s = 1/t_s; %采樣頻率
f_x = 0:f_s/(Sampling_points -1):f_s;  %注意這里和橫坐標(biāo)頻率對(duì)應(yīng)上了，頻率分辨率就是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頻率分量移到坐標(biāo)中心
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結(jié)果的前半部分

奈奎斯特定理的作用

奈奎斯特定理在頻域分析中扮演著核心角色。它指出，信號(hào)的采樣速度必須高于其最高頻率的兩倍。比如，對(duì)于雷達(dá)信號(hào)，如果雷達(dá)發(fā)射的最高頻率是10MHz，那么采樣速度至少要是20MHz。如果低于這個(gè)標(biāo)準(zhǔn)，就會(huì)導(dǎo)致頻譜重疊。實(shí)驗(yàn)表明，如果采樣速度不夠，采集到的頻譜圖中會(huì)出現(xiàn)錯(cuò)誤頻率，這會(huì)影響到所有頻域處理結(jié)果的精確度。

換個(gè)角度審視，此定理明確了頻域分析中基礎(chǔ)采樣準(zhǔn)則的設(shè)定。比如，在處理地震波信號(hào)時(shí)，若不依照奈奎斯特定理，低頻地震信號(hào)可能會(huì)與高頻噪聲混淆，進(jìn)而造成地層結(jié)構(gòu)關(guān)鍵信息的識(shí)別困難。

功率譜的求法

計(jì)算功率譜有倆種方法。一種是用傅立葉變換的平方除以區(qū)間長(zhǎng)度來算，這方法挺簡(jiǎn)單，容易理解。在MATLAB編程里，我們根據(jù)這個(gè)公式就能得到結(jié)果。例如，針對(duì)描述粒子運(yùn)動(dòng)軌跡的離散信號(hào)，用這個(gè)方法算功率譜，就能看出粒子在不同頻率上的能量分布情況。

還有一法是采用自相關(guān)函數(shù)執(zhí)行傅里葉轉(zhuǎn)換。從理論層面分析，這兩種手段的成果本應(yīng)相同。但實(shí)際操作中，例如在分析含噪電子設(shè)備輸出信號(hào)時(shí)，我們發(fā)現(xiàn)自相關(guān)函數(shù)法在降噪效果上更為突出，繪制的曲線也更加平滑。這好比在喧鬧場(chǎng)合努力辨識(shí)對(duì)方說話，運(yùn)用自相關(guān)函數(shù)傅里葉轉(zhuǎn)換這一技術(shù)，能更高效地去除噪音干擾，進(jìn)而更清晰地呈現(xiàn)出信號(hào)的功率譜特征。

Fs = 1000;
nfft = 1000;  %fft采樣點(diǎn)數(shù)

%產(chǎn)生序列
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('加噪信號(hào)');xlim([0 1000]);grid on
%FFT
Y = 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
%自相關(guān)結(jié)果
cxn = xcorr(xn,'unbiased');  %計(jì)算自相關(guān)函數(shù)
%自相關(guān)法
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

直接法與周期圖法對(duì)比

編寫代碼時(shí)，我們運(yùn)用直接法，通過傅立葉變換的平方除以區(qū)間長(zhǎng)度來求得功率譜。這一方法與MATLAB中的periodogram函數(shù)所得結(jié)果相同。在分析機(jī)械振動(dòng)信號(hào)時(shí)，我們對(duì)比了這兩種計(jì)算功率譜的方法。結(jié)果顯示，在數(shù)據(jù)量不多且信號(hào)頻率成分簡(jiǎn)單的情況下，兩種方法得出的結(jié)果非常接近。這一發(fā)現(xiàn)為我們?cè)诰帉戭l域分析代碼時(shí)提供了更多選擇，使我們能根據(jù)實(shí)際情況挑選出更易理解或計(jì)算速度更快的方案。

工業(yè)監(jiān)測(cè)系統(tǒng)若要依賴實(shí)時(shí)功率譜來預(yù)測(cè)故障，這一點(diǎn)極為關(guān)鍵。以風(fēng)力發(fā)電機(jī)的振動(dòng)監(jiān)控為例，我們需根據(jù)硬件條件和信號(hào)特點(diǎn)，巧妙選擇計(jì)算功率譜的方法。最重要的是，要能快速且精確地識(shí)別出潛在的危險(xiǎn)頻率信號(hào)。

倒頻譜的計(jì)算與含義

實(shí)倒頻譜函數(shù)用于計(jì)算倒頻譜，MATLAB手冊(cè)中指出其計(jì)算步驟是先將信號(hào)轉(zhuǎn)換為頻譜，接著轉(zhuǎn)為對(duì)數(shù)形式，最后執(zhí)行傅立葉逆變換。但倒頻譜的定義是將信號(hào)轉(zhuǎn)成功率譜，隨后進(jìn)行對(duì)數(shù)處理，再進(jìn)行傅立葉逆變換。兩者間的差異可能源于功率譜是頻譜值的平方，對(duì)數(shù)處理后，平方變成了系數(shù)2，但這對(duì)于后續(xù)計(jì)算影響微乎其微，所以近似結(jié)果并無差異。

在仿真實(shí)驗(yàn)里，我們需要留意倒頻譜的作用，需要親手制作一組調(diào)頻信號(hào)。舉例來說，對(duì)那些高頻（如50Hz、100Hz、200Hz）和低頻（如5Hz、10Hz、20Hz）的信號(hào)進(jìn)行調(diào)整，接著分別畫出低頻、高頻和調(diào)頻信號(hào)的時(shí)域圖和頻譜圖。這種調(diào)頻信號(hào)的處理方法，在通信信號(hào)分析領(lǐng)域應(yīng)用很廣，有助于提高信號(hào)傳輸?shù)男屎头€(wěn)定性。

倒頻譜在信號(hào)處理中的體現(xiàn)

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);
end
subplot(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/sf:(nfft-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('時(shí)間(s)');
ylabel('幅值');
grid on;
subplot(3,3,8)
yy = real(ifft(log(abs(fft(y))))); %信號(hào)→傅里葉→對(duì)數(shù)→傅里葉逆變換
plot(t(nn),abs(yy(nn)));
title('real(ifft(log(abs(fft(y)))))');ylim([0 0.3]);
xlabel('時(shí)間(s)');
ylabel('幅值');
grid on;

以特定調(diào)制信號(hào)為例，例如某些信號(hào)中帶有20Hz、10Hz和5Hz等低頻部分。在常規(guī)的FFT_y分析中，這些低頻部分以邊緣頻帶的形式出現(xiàn)，難以準(zhǔn)確判斷其頻率。但在倒頻譜分析中，這些部分卻很容易被識(shí)別。這好比在黑夜中，一雙隱藏的眼睛在FFT_y分析中難以發(fā)現(xiàn)，而在倒頻譜分析中卻能立刻被“照亮”，顯示出其特征。在無線電通信的信號(hào)分析中，倒頻譜分析這種處理類似問題的能力對(duì)工程師來說非常寶貴，它能夠幫助他們快速發(fā)現(xiàn)并解讀微弱的低頻信號(hào)。

關(guān)于項(xiàng)目操作中選用頻域分析方法的問題，我想聽聽大家的意見：面對(duì)不同項(xiàng)目的具體特點(diǎn)，你們是如何做出選擇的？歡迎在評(píng)論區(qū)留下你們的想法。若這篇文章對(duì)您有所啟發(fā)，不妨點(diǎn)贊或分享給更多人。