Matlab利用離散傅里葉變換DFT進(jìn)行頻譜分析的步驟

信號在頻域能夠呈現(xiàn)出時域不易發(fā)現(xiàn)的性質(zhì)和規(guī)律，傅里葉變換是將信號從時域變換到頻域，便于在頻域?qū)π盘柕奶匦赃M(jìn)行分析。離散傅里葉變換 (DFT)，是傅里葉變換在時域和頻域上的離散呈現(xiàn)形式，通俗的說就是將經(jīng)過采樣的有限長度時域離散采樣序列變換為等長度的頻域離散采樣序列，通過對變換得到的頻域采樣序列進(jìn)行適當(dāng)?shù)膿Q算和處理，可以得到信號的頻譜（頻率-幅值曲線和頻率-相位曲線）。

1. 定義

離散傅里葉變換 (DFT)的定義為：

式中，x(n)為時域離散采樣序列（通常為實(shí)數(shù)序列），N為時域離散采樣序列x(n)的長度，X(k)為頻域離散采樣序列（通常為復(fù)數(shù)序列）。

快速傅里葉變換（FFT）是離散傅里葉變換（DFT）的一種快速算法，F(xiàn)FT的計算結(jié)果與DFT完全相同，但FFT相對于DFT減小了計算量、節(jié)約計算資源消耗，能夠適應(yīng)在線計算，因此實(shí)際DFT都是通過FFT算法來求得結(jié)果。

Matlab軟件自帶fft函數(shù)實(shí)現(xiàn)快速傅里變換算法，但是光使用fft并不能直接得到信號的頻譜，還需要解決以下問題：

• 幅值變換 ：X(k)序列的幅值大小與參與變換的時域序列x(n)長度N有關(guān)，變換后的幅值|X(k)|需要乘以2/N得到真實(shí)幅值；
• 有效頻率區(qū)域 ：X(k)序列由兩部分共軛復(fù)數(shù)序列組成（復(fù)數(shù)共軛表示幅值相等、相位相反），相當(dāng)于只有一半的復(fù)數(shù)序列是獨(dú)立有效的，這部分復(fù)數(shù)序列對應(yīng)0~fs/2的頻率區(qū)域（fs為時域離散采樣序列x(n)的采樣頻率）。
• ** 直流信號的處理** ：直流信號幅值（對應(yīng)頻率0Hz）為兩部分共軛復(fù)數(shù)序列在頻率0Hz處的加和，其真實(shí)幅值再乘以2/N后還需要再除以2得到真實(shí)的直流信號幅值。

2. 函數(shù)

初學(xué)的朋友若不理解上述變換和處理技巧，很難得到正確的頻譜圖。為此作者在fft函數(shù)的基礎(chǔ)上，使用Matlab開發(fā)了函數(shù) DFT.m ，通過函數(shù)來實(shí)現(xiàn)上述幅值變換、有效頻率區(qū)域和直流信號的處理，能夠直接分析出給定離散信號x(n)的幅值譜和相位譜，函數(shù)簡單、易用、通用性好。

function [f,X_m,X_phi] = DFT(xn,ts,N,drawflag)
% [f,X_m,X_phi] = DFT(xn,ts,N,drawflag) 離散序列的快速傅里葉變換，時域轉(zhuǎn)換為頻域
% 輸入  xn為離散序列 為向量  
%       ts為序列的采樣時間/s
%       N為FFT變換的點(diǎn)數(shù)，默認(rèn)為xn的長度  
%       drawflag為繪圖標(biāo)識位，取0時不繪圖，其余非0值時繪圖，默認(rèn)為繪圖
% 輸出 f為頻率向量
%      X_m為幅值向量
%      X_phi為相位向量，單位為°
% 注意計算出來的0頻分量(直流分量應(yīng)該除以2)  直流分量的符號應(yīng)結(jié)合相位圖來確定

3. 演示

3.1 單頻正弦信號（整數(shù)周期采樣）

%% Eg 1 單頻正弦信號
ts = 0.01;
t = 0:ts:1;
A = 1.5;       % 幅值  
f = 2;         % 頻率
w = 2*pi*f;    % 角頻率
phi = pi/3;    % 初始相位 
x = A*cos(w*t+phi);   % 時域信號
figure
plot(t,x)
xlabel('時間/s')
ylabel('時域信號x(t)')
% DFT變換將時域轉(zhuǎn)換到頻域,并繪制頻譜圖
[f,X_m,X_phi] = DFT(x,ts);

結(jié)果 ：

正弦信號頻率為2Hz，頻譜分析頻率為1.98Hz
正弦信號幅值為1.5，頻譜分析幅值為1.495
正弦信號相位為60°，頻譜分析相位為63.32°

**3.2 **單頻正弦信號（非整數(shù)周期采樣）

%% Eg 2 單頻正弦信號(非整數(shù)周期采樣)
ts = 0.01;
t = 0:ts:1;
A = 1.5;       % 幅值  
f = 1.5;         % 頻率
w = 2*pi*f;    % 角頻率
phi = pi/3;    % 初始相位 
x = A*cos(w*t+phi);   % 時域信號
figure
plot(t,x)
xlabel('時間/s')
ylabel('時域信號x(t)')
% DFT變換將時域轉(zhuǎn)換到頻域,并繪制頻譜圖
[f,X_m,X_phi] = DFT(x,ts);

結(jié)果 ：

正弦信號頻率為1.5Hz，頻譜分析頻率為0.99Hz、1.98Hz
正弦信號幅值為1.5，頻譜分析幅值為1.034、0.923
正弦信號相位為60°，頻譜分析相位為160.93°、-33.76°

總結(jié)：

DFT變換后頻率序列的最小單位刻度為f s /N（此例為1Hz），非整數(shù)周期采樣時關(guān)心信號的頻率（此例為1.5Hz）不是頻率分辨率f s /N的正整數(shù)倍，那這個頻率成分信號會由前后兩個正整數(shù)倍的頻率成分信號（此例為1Hz和2Hz）的線性組合來替代，這就是頻譜泄漏現(xiàn)象，非周期采樣時某頻率成分信號向兩側(cè)頻率分辨率正整數(shù)倍的頻點(diǎn)泄漏。實(shí)際頻譜分析時并不清楚所關(guān)心的頻率點(diǎn)精確值，避免此問題的一個解決方法是，取更多的點(diǎn)參加DFT，即時域序列x(n)長度N值取長一些，讓頻率分辨率f s /N很小，以減小頻譜泄漏現(xiàn)象。再看下例。

%% Eg 2 單頻正弦信號(非整數(shù)周期采樣)
ts = 0.01;
t = 0:ts:1;
A = 1.5;       % 幅值  
f = 1.5;         % 頻率
w = 2*pi*f;    % 角頻率
phi = pi/3;    % 初始相位 
x = A*cos(w*t+phi);   % 時域信號
figure
plot(t,x)
xlabel('時間/s')
ylabel('時域信號x(t)')
% DFT變換將時域轉(zhuǎn)換到頻域,并繪制頻譜圖
[f,X_m,X_phi] = DFT(x,ts);

結(jié)果： 頻譜泄漏情況大為改善，采樣點(diǎn)繼續(xù)增多時，頻譜泄漏會進(jìn)一步減小。

正弦信號頻率為1.5Hz，頻譜分析頻率主要成分為1.46Hz、
正弦信號幅值為1.5，頻譜分析頻率主要成分對應(yīng)幅值為1.41
正弦信號相位為60°，頻譜分析頻率主要成分對應(yīng)相位為89.5°

3.3 含有直流分量的單頻正弦信號

%% Eg 3 含有直流分量的單頻正弦信號
ts = 0.01;
t = 0:ts:1;
A = 1.5;       % 幅值  
f = 5;         % 頻率
w = 2*pi*f;    % 角頻率
phi = pi/6;    % 初始相位 
x = 0.5 + A*cos(w*t+phi);   % 時域信號,帶有直流偏移0.5
figure
plot(t,x)
xlabel('時間/s')
ylabel('時域信號x(t)')
% DFT變換將時域轉(zhuǎn)換到頻域,并繪制頻譜圖
[f,X_m,X_phi] = DFT(x,ts);

結(jié)果 ：

正弦信號頻率為5Hz，頻譜分析頻率為4.95Hz
正弦信號幅值為1.5，頻譜分析幅值為1.498
正弦信號相位為30°，頻譜分析相位為38.66°
正弦信號直流分量0.5，頻譜分析直流分量為0.51

**3.4 **正弦復(fù)合信號

%% Eg 4 正弦復(fù)合信號
ts = 0.01;
t = 0:ts:2;
A = [1.5 1 0.5 0.2];    % 幅值  
f = [3 6 9 15];         % 頻率
w = 2*pi*f;             % 角頻率
phi = (1:4)*pi/4;       % 初始相位 
x = -0.5 + A(1)*cos(w(1)*t+phi(1)) + A(2)*cos(w(2)*t+phi(2)) + A(3)*cos(w(3)*t+phi(3)) + A(4)*cos(w(4)*t+phi(4));     % 時域信號
figure
plot(t,x)
xlabel('時間/s')
ylabel('時域信號x(t)')
% DFT變換將時域轉(zhuǎn)換到頻域,并繪制頻譜圖
[f,X_m,X_phi] = DFT(x,ts);

結(jié)果 ：

正弦信號頻率為3、6、9、15Hz，頻譜分析頻率為2.985、5.97、8.96、14.93Hz
正弦信號幅值為1.5、1、0.5、0.2，頻譜分析幅值為1.499、0.989、0.485、0.192
正弦信號相位為45°、90°、135°、180°，頻譜分析相位為50.6°、101.5°、152.9°、210°
正弦信號直流分量-0.5，頻譜分析直流分量為-0.497

注意 ：頻率為0Hz時對應(yīng)的直流信號的幅值的正負(fù)號，是通過零頻相位來確定的，相位為0°表示幅值為正，相位為180°表示幅值為負(fù)。

**3.5 **含有隨機(jī)干擾的正弦信號

%% Eg 5 含有隨機(jī)干擾的正弦信號
ts = 0.01;
t = 0:ts:2;
A = [1 0.5];    % 幅值  
f = [3 10];         % 頻率
w = 2*pi*f;             % 角頻率
phi = (1:2)*pi/3;       % 初始相位 
x =  A(1)*cos(w(1)*t+phi(1)) + A(2)*cos(w(2)*t+phi(2)) + 0.8*(rand(size(t))-0.5);     % 時域信號
figure
plot(t,x)
xlabel('時間/s')
ylabel('時域信號x(t)')
% DFT變換將時域轉(zhuǎn)換到頻域,并繪制頻譜圖
[f,X_m,X_phi] = DFT(x,ts);

結(jié)果 ：

正弦信號頻率為3、10Hz，頻譜分析頻率為2.985、9.95Hz
正弦信號幅值為1、0.5，頻譜分析幅值為0.978、0.456
正弦信號相位為60°、135°，頻譜分析相位為65.1°、139.8°

**3.6 **實(shí)際案例

load data
ts = 0.001;
x = Jsd;
t = [0:length(x)-1]*ts;
figure
plot(t,x)
xlabel('時間/s')
ylabel('時域信號x(t)')
% DFT變換將時域轉(zhuǎn)換到頻域,并繪制頻譜圖
[f,X_m,X_phi] = DFT(x,ts);

結(jié)果 ：

頻譜分析主要頻率成分為18.996、37.992Hz
頻譜分析主要頻率成分對應(yīng)幅值為1.741、1.117

該項(xiàng)目為作者在強(qiáng)振環(huán)境下測得加速度信號，加速度是機(jī)械結(jié)構(gòu)周期運(yùn)動激勵產(chǎn)生，需要通過頻譜分析獲取機(jī)械結(jié)構(gòu)周期運(yùn)動的頻率。由于噪聲幅度遠(yuǎn)大于有效信號幅度，信號的信噪比很低，從時域上很難辨別機(jī)械結(jié)構(gòu)周期運(yùn)動的頻率。但經(jīng)過DFT后，從頻域上可以看出信號的主要頻率成分為19Hz和其倍頻38Hz，可以判斷機(jī)械結(jié)構(gòu)周期運(yùn)動的頻率為19Hz，38Hz為結(jié)構(gòu)響應(yīng)的非線性特性所產(chǎn)生的倍頻。