使用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');

運行結(jié)果如下所示：

image

快速傅里葉

matlab中的快速傅里葉有兩種調(diào)用形式：

y=fft(x)吗氏。x是原始信號乌昔，y是變換之后的信號。y與x有相同的長度尾组。
y=fft(x,N)忙芒。x，y定義如上讳侨，N是正整數(shù)呵萨，表示進行N點快速傅里葉變換。如果x長度小于N爷耀，則對x補零甘桑，使之與N相等；反則，則對x進行截取跑杭。

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

一般而言德谅，N點fft的結(jié)果y爹橱，在 $n=N/2+1$ 處對應(yīng)的頻率為最高采樣率的一半，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;%橫軸各點對應(yīng)的頻率值

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

結(jié)果圖省略椭盏。

濾波

利用快速傅里葉簡單濾波

下例是將振幅為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;%橫軸各點對應(yīng)的頻率值

%濾波
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可能是復(fù)數(shù)掏颊，對其取實部得到信號
tt = [0:length(xx)-1]/fs;%橫軸各點對應(yīng)的時間

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

結(jié)果如下圖

image

幾個簡單的函數(shù)

xi=filtic(B,A,ys)糟红。B、A分別為系統(tǒng)z變換后的傳遞函數(shù)的分子和分母多項式的系數(shù)向量乌叶，ys是系統(tǒng)的初始輸出狀態(tài)盆偿，xi為對應(yīng)的初始條件下輸入序列。
yn0=filter(B,A,xn)准浴。B事扭、A定義同上，xn是系統(tǒng)的輸入信號乐横，yn0為系統(tǒng)的零狀態(tài)響應(yīng)求橄。
yn=filter(B,A,xn,xi)。B晰奖、A谈撒、xn、xi定義同上匾南，yn為系統(tǒng)全響應(yīng)啃匿。

模擬濾波器

以巴特沃斯低通濾波器為例，說明調(diào)用方法蛆楞。

[btt1,ctt1] = butter(N,wn,'s');%1.調(diào)用函數(shù)生成濾波器系數(shù)
H = [tf(btt1,ctt1)];%濾波器的傳遞函數(shù)
t = (0:n-1)./fs;%時域信號橫軸的坐標溯乒，n為長度，fs為采樣率
s1 = lsim(H,a1,t);%2.濾波

說明：

[btt1,ctt1] = butter(N,wn,'s');豹爹。N是濾波器的階數(shù)裆悄，wn是截止頻率（是弧度值，如果截止頻率要求為500Hz臂聋，則 $wn=500*2*\pi$ ）光稼』蚰希可以直接給定，亦可以根據(jù)參數(shù)由buttord`函數(shù)計算得到艾君。's'表示模擬濾波器采够。btt1、ctt1分別表示濾波器在拉普拉斯域中傳遞函數(shù)的分子冰垄、分母多項式的系數(shù)蹬癌。
s1 = lsim(H,a1,t)。H是模擬濾波器的傳遞函數(shù)虹茶，a1表示待濾波信號逝薪，t是信號的橫坐標，s1是濾波后的信號蝴罪。

其他說明：

這里僅以低通濾波器為例董济，其他巴特沃斯濾波器如高通、帶通洲炊、帶阻調(diào)用方式類似感局，只是函數(shù)butter的參數(shù)略有不同，請參看matlab關(guān)于butter函數(shù)的介紹暂衡。（在matlab中執(zhí)行help butter)
其他濾波器，如橢圓濾波器等崖瞭，使用方式類似狂巢，只是函數(shù)名稱不同。

數(shù)字濾波器

以巴特沃斯低通濾波器為例书聚，說明調(diào)用方法唧领。

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

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

說明:

[btt,ctt] = butter(N,wn)。N是濾波器階數(shù)雌续，wn是相對截止頻率斩个，比如最高采樣率為Fs，要求的截止頻率為fs驯杜，則 $wn=fs/Fs$ 受啥。可以直接給定鸽心，亦可以根據(jù)參數(shù)由buttord函數(shù)計算得到滚局。注意，這里沒有參數(shù)'s'顽频。btt藤肢、ctt分別表示濾波器在z域中傳遞函數(shù)的分子、分母多項式的系數(shù)糯景。
Signal_Filter=filter(btt,ctt,a1)嘁圈。btt省骂、ctt與之前定義相同，a1是待濾波信號最住，Signal_Filter是濾波之后的信號钞澳。
[btt1,ctt1]=bilinear(btt1,ctt1,Fs)。是使用雙線性法將模擬濾波器在拉普拉斯域中的系數(shù)轉(zhuǎn)換成數(shù)字濾波器在z域中的系數(shù)温学，F(xiàn)s是采樣率略贮。

其他說明：

這里僅以低通濾波器為例，其他巴特沃斯濾波器如高通仗岖、帶通逃延、帶阻調(diào)用方式類似，只是函數(shù)butter的參數(shù)略有不同轧拄，請參看matlab關(guān)于butter函數(shù)的介紹揽祥。（在matlab中執(zhí)行help butter)
其他濾波器，如橢圓濾波器等檩电，使用方式類似拄丰，只是函數(shù)名稱不同。