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

1. 定義

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

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

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

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

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

2. 函數(shù)

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

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

3. 演示

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

%% Eg 1 單頻正弦信號(hào)
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);   % 時(shí)域信號(hào)
figure
plot(t,x)
xlabel('時(shí)間/s')
ylabel('時(shí)域信號(hào)x(t)')
% DFT變換將時(shí)域轉(zhuǎn)換到頻域,并繪制頻譜圖
[f,X_m,X_phi] = DFT(x,ts);

結(jié)果 ：

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

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

%% Eg 2 單頻正弦信號(hào)(非整數(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);   % 時(shí)域信號(hào)
figure
plot(t,x)
xlabel('時(shí)間/s')
ylabel('時(shí)域信號(hào)x(t)')
% DFT變換將時(shí)域轉(zhuǎn)換到頻域,并繪制頻譜圖
[f,X_m,X_phi] = DFT(x,ts);

結(jié)果 ：

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

總結(jié)：

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

%% Eg 2 單頻正弦信號(hào)(非整數(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);   % 時(shí)域信號(hào)
figure
plot(t,x)
xlabel('時(shí)間/s')
ylabel('時(shí)域信號(hào)x(t)')
% DFT變換將時(shí)域轉(zhuǎn)換到頻域,并繪制頻譜圖
[f,X_m,X_phi] = DFT(x,ts);

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

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

3.3 含有直流分量的單頻正弦信號(hào)

%% Eg 3 含有直流分量的單頻正弦信號(hào)
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);   % 時(shí)域信號(hào),帶有直流偏移0.5
figure
plot(t,x)
xlabel('時(shí)間/s')
ylabel('時(shí)域信號(hào)x(t)')
% DFT變換將時(shí)域轉(zhuǎn)換到頻域,并繪制頻譜圖
[f,X_m,X_phi] = DFT(x,ts);

結(jié)果 ：

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

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

%% Eg 4 正弦復(fù)合信號(hào)
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));     % 時(shí)域信號(hào)
figure
plot(t,x)
xlabel('時(shí)間/s')
ylabel('時(shí)域信號(hào)x(t)')
% DFT變換將時(shí)域轉(zhuǎn)換到頻域,并繪制頻譜圖
[f,X_m,X_phi] = DFT(x,ts);

結(jié)果 ：

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

注意 ：頻率為0Hz時(shí)對(duì)應(yīng)的直流信號(hào)的幅值的正負(fù)號(hào)，是通過(guò)零頻相位來(lái)確定的，相位為0°表示幅值為正，相位為180°表示幅值為負(fù)。

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

%% Eg 5 含有隨機(jī)干擾的正弦信號(hào)
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);     % 時(shí)域信號(hào)
figure
plot(t,x)
xlabel('時(shí)間/s')
ylabel('時(shí)域信號(hào)x(t)')
% DFT變換將時(shí)域轉(zhuǎn)換到頻域,并繪制頻譜圖
[f,X_m,X_phi] = DFT(x,ts);

結(jié)果 ：

正弦信號(hào)頻率為3、10Hz，頻譜分析頻率為2.985、9.95Hz
正弦信號(hào)幅值為1、0.5，頻譜分析幅值為0.978、0.456
正弦信號(hào)相位為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('時(shí)間/s')
ylabel('時(shí)域信號(hào)x(t)')
% DFT變換將時(shí)域轉(zhuǎn)換到頻域,并繪制頻譜圖
[f,X_m,X_phi] = DFT(x,ts);

結(jié)果 ：

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

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