傅立葉分析

公式法

下例 是將振幅為1的5Hz正弦波和振幅為0.5的10Hz正弦波相加之後進行傅立葉分析。

clear all
N=512;
dt=0.02;
n=0:N-1;
t=n*dt;
x=sin(2*pi*5*t)+0.5*sin(2*pi*10*t);%生成和訊號

%傅立葉變換
m = floor(N/2)+1;
a=zeros(1,m);
b=zeros(1,m);

for k=0:m-1
    for ii=0:N-1
        a(k+1) = a(k+1)+2/N*x(ii+1)*cos(2*pi*k*ii/N);
        b(k+1) = b(k+1)+2/N*x(ii+1)*sin(2*pi*k*ii/N);
    end
    c(k+1)=sqrt(a(k+1).^2+b(k+1).^2);
end

%傅立葉逆變換
if(mod(N,2) ~=1)
    a(m)=a(m)/2;
end
for ii=0:N-1
    xx(ii+1)=a(1)/2;
    for k=1:m-1
        xx(ii+1)=xx(ii+1)+a(k+1)*cos(2*pi*k*ii/N)+b(k+1)*sin(2*pi*k*ii/N);
    end
end

%繪圖
subplot(3,1,1),plot(t,x,'LineWidth',2);title('原始訊號'),xlabel('時間/s');
subplot(3,1,2),plot((0:m-1)/(N*dt),c,'LineWidth',2);title('傅立葉變換'),xlabel('頻率/Hz');
subplot(3,1,3),plot((0:N-1)*dt,xx,'LineWidth',2);title('合成訊號'),xlabel('時間/s');
 

執行結果如下所示：

快速傅立葉

matlab中的快速傅立葉有兩種呼叫形式：

• y=fft(x)。x是原始訊號，y是變換之後的訊號。y與x有相同的長度。
• y=fft(x,N)。x，y定義如上，N是正整數，表示進行N點快速傅立葉變換。如果x長度小於N，則對x補零，使之與N相等；反則，則對x進行擷取。

對應的逆變換有兩種，分別為x=ifft(y)和x=ifft(y.N)。

一般而言，N點fft的結果y，在$n=N/2+1$處對應的頻率為最高取樣率的一半，y的後一半與前一半對稱。

下例 是將振幅為1的5Hz正弦波和振幅為0.5的10Hz正弦波相加之後進行傅立葉分析。

clc;clear;
fs=30;
N=512;
n=0:N-1;
t=n/fs;
x=sin(2*pi*5*t)+0.5*sin(2*pi*10*t);

%快速傅立葉
y=fft(x,N);
mag = abs(y);%振幅
ang = angle(y);%相位
f = n*fs/N;%橫軸各點對應的頻率值

%逆變換
xx = ifft(y);
xx = real(xx);%計算誤差使得xx可能是複數，對其取實部得到訊號
tt = [0:length(xx)-1]/fs;%橫軸各點對應的時間
 

結果圖省略。

濾波

利用快速傅立葉簡單濾波

下例是將振幅為1的5Hz正弦波和振幅為0.5的10Hz正弦波相加之後，濾除8Hz以上的訊號。

clc;clear;
fs=30;%取樣率
N=256;
n=0:N-1;
t=n/fs;
dt=1/fs;
x=sin(2*pi*5*t)+0.5*sin(2*pi*10*t);

%快速傅立葉
y=fft(x,N);
mag = abs(y)*2/N;%振幅
ang = angle(y);%相位
f = n*fs/N;%橫軸各點對應的頻率值

%濾波
nn=length(y);
yy = zeros(1,nn);
for m=0:nn-1
    if(m/(nn*dt)>8)&(m/(N*dt)<(1/dt-8))
        yy(m+1)=0;
    else
        yy(m+1)=y(m+1);
    end
end

%逆變換
xx = ifft(yy);
xx = real(xx);%計算誤差使得xx可能是複數，對其取實部得到訊號
tt = [0:length(xx)-1]/fs;%橫軸各點對應的時間

%繪圖
subplot(2,1,1),plot(t,x,'LineWidth',2);title('原始訊號'),xlabel('時間/s');
subplot(2,1,2),plot(tt,xx,'LineWidth',2);title('濾波後的訊號'),xlabel('時間/s');

 

結果如下圖

幾個簡單的函式

• xi=filtic(B,A,ys)。B、A分別為系統z變換後的傳遞函式的分子和分母多項式的係數向量，ys是系統的初始輸出狀態，xi為對應的初始條件下輸入序列。
• yn0=filter(B,A,xn)。B、A定義同上，xn是系統的輸入訊號，yn0為系統的零狀態響應。
• yn=filter(B,A,xn,xi)。B、A、xn、xi定義同上，yn為系統全響應。

模擬濾波器

以巴特沃斯低通濾波器為例，說明呼叫方法。

[btt1,ctt1] = butter(N,wn,'s');%1.呼叫函式生成濾波器係數
H = [tf(btt1,ctt1)];%濾波器的傳遞函式
t = (0:n-1)./fs;%時域訊號橫軸的座標，n為長度，fs為取樣率
s1 = lsim(H,a1,t);%2.濾波

說明：

1. [btt1,ctt1] = butter(N,wn,'s');。N是濾波器的階數，wn是截止頻率（是弧度值，如果截止頻率要求為500Hz，則$wn=5002\pi$）。可以直接給定，亦可以根據引數由buttord`函式計算得到。's'表示模擬濾波器。btt1、ctt1分別表示濾波器在拉普拉斯域中傳遞函式的分子、分母多項式的係數。
2. s1 = lsim(H,a1,t)。H是模擬濾波器的傳遞函式，a1表示待濾波訊號，t是訊號的橫座標，s1是濾波後的訊號。

其他說明：

• 這裡僅以低通濾波器為例，其他巴特沃斯濾波器如高通、帶通、帶阻呼叫方式類似，只是函式butter的引數略有不同，請參看matlab關於butter函式的介紹。（在matlab中執行help butter)
• 其他濾波器，如橢圓濾波器等，使用方式類似，只是函式名稱不同。

數字濾波器

以巴特沃斯低通濾波器為例，說明呼叫方法。

%方式一：直接設計
[btt,ctt] = butter(N,wn);%1.生成數字濾波器
Signal_Filter=filter(btt,ctt,a1);%2.濾波

%方式二：模擬濾波器轉數字濾波器
[btt1,ctt1] = butter(N,Wn,'s');
[btt1,ctt1]=bilinear(btt1,ctt1,Fs);%3.模擬濾波器轉數字濾波器
Signal_Filter=filter(btt,ctt,a1);

說明:

1. [btt,ctt] = butter(N,wn)。N是濾波器階數，wn是相對截止頻率，比如最高取樣率為Fs，要求的截止頻率為fs，則$wn=fs/Fs$ 。可以直接給定，亦可以根據引數由buttord函式計算得到。注意，這裡沒有引數's'。btt、ctt分別表示濾波器在z域中傳遞函式的分子、分母多項式的係數。
2. Signal_Filter=filter(btt,ctt,a1)。btt、ctt與之前定義相同，a1是待濾波訊號，Signal_Filter是濾波之後的訊號。
3. [btt1,ctt1]=bilinear(btt1,ctt1,Fs)。是使用雙線性法將模擬濾波器在拉普拉斯域中的係數轉換成數字濾波器在z域中的係數，Fs是取樣率。

其他說明：

• 這裡僅以低通濾波器為例，其他巴特沃斯濾波器如高通、帶通、帶阻呼叫方式類似，只是函式butter的引數略有不同，請參看matlab關於butter函式的介紹。（在matlab中執行help butter)
• 其他濾波器，如橢圓濾波器等，使用方式類似，只是函式名稱不同。