使用matlab進行傅里葉分析和濾波

傅里葉分析

公式法

下例 是將振幅為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)
• 其他濾波器，如橢圓濾波器等，使用方式類似，只是函數名稱不同。