MATLAB訊號處理——數字濾波器的設計

有覺悟的小韭菜發表於2020-10-13

一、模擬FIR濾波器的設計函式(濾波器的階數是指在濾波器的傳遞函式中有幾個極點)

1.模擬濾波器的階數選擇函式(以巴特沃斯模擬濾波器為例)

(1)函式呼叫:

[N,Wn]=buttord(wp,ws,Rp,Rs,'s')

(2)例子:

wp=[0.2*pi 0.3*pi];              % 設定通帶頻率
ws=[0.1*pi 0.4*pi];              % 設定阻帶頻率
Rp=1; Rs=20;                     % 設定波紋係數
% 巴特沃斯濾波器設計
[N,Wn]=buttord(wp,ws,Rp,Rs,'s'); % 求巴特沃斯濾波器階數

其中N是濾波器最小階數,Wn是等效低通濾波器的截止頻率

2.在MATLAB中設計模擬濾波器的函式

(1)函式呼叫

[bb,ab]=butter(N,Wn,'s');        % 求巴特沃斯濾波器係數

(2)例子

[bb,ab]=butter(N,Wn,'s');        % 求巴特沃斯濾波器係數

3.濾波器頻率分析函式

(1)freqs函式

①函式呼叫:
在這裡插入圖片描述
②例子:

clear all; close all; clc;
wp=[0.2*pi 0.3*pi];              % 設定通帶頻率
ws=[0.1*pi 0.4*pi];              % 設定阻帶頻率
Rp=1; Rs=20;                     % 設定波紋係數
% 巴特沃斯濾波器設計
[N,Wn]=buttord(wp,ws,Rp,Rs,'s'); % 求巴特沃斯濾波器階數
fprintf('巴特沃斯濾波器 N=%4d\n',N) % 顯示濾波器階數
[bb,ab]=butter(N,Wn,'s');        % 求巴特沃斯濾波器係數
W=0:0.01:2;                      % 設定模擬頻率
[Hb,wb]=freqs(bb,ab,W);          % 求巴特沃斯濾波器頻率響應
plot(wb/pi,20*log10(abs(Hb)),'b')% 作圖

在這裡插入圖片描述

4.模擬濾波器的數字化方法

在這裡插入圖片描述

二、IIR數字濾波器的設計函式(本文Wp為模擬角頻率,wp為數字角頻率)

1.IIR數字濾波器階數選擇函式

(1)函式呼叫:

[N,Wn]=buttord(wp,ws,Rp,Rs)

2.在MATLAB中確定濾波器係數

(1)函式呼叫

[bb,ab]=butter(N,Wn,'s');        % 求巴特沃斯濾波器係數

3.模擬濾波器數字化

(1)脈衝響應不變法

①具體步驟:在這裡插入圖片描述
②函式呼叫:
在這裡插入圖片描述
③脈衝響應不變法造成的混疊失真
在這裡插入圖片描述
④脈衝響應不變法的最大特點是頻率響應的混疊響應,所以脈衝響應不變法只適用於衰減特性很好的低通或帶通濾波器,在高頻衰減越快,混疊效應越小。
⑤使用impinvar函式進行脈衝響應不變法的例項:

clear all; clc; close all;
bs=[1,1];as=[1,5,6];            % 系統分子分母系數向量
Fs=10; T=1/Fs;                  % 取樣頻率和取樣間隔
% 呼叫impinvar函式計算數字濾波器係數
[Bd,Ad]=impinvar(bs,as,Fs);

⑥不能使用impinvar函式進行脈衝響應不變法的原因以及解決方法:
1°濾波器型別選擇不合適
2°過渡帶選擇太窄。
3°選擇由模擬濾波器轉換成數字濾波器的方法不合適

(2)雙線性Z變換法

①步驟
1°根據要求進行頻率預畸處理
在這裡插入圖片描述
2°根據預畸變後的頻率求濾波器係數
②例:

clear all; close all; clc;
fp=100; fs=200;                  % 設定通帶和阻帶
Fs=1000;                         % 取樣頻率
Rp=2; Rs=40;                     % 通帶波紋和阻帶衰減
wp=fp*2*pi/Fs;                   % 把通帶和阻帶設為角頻率
ws=fs*2*pi/Fs;
T=1;                             % T=1
Ts=1/Fs;                         % Ts=1/Fs
Wp=2/Ts*tan(wp/2);               % 把通帶和阻帶按Fs進行預畸
Ws=2/Ts*tan(ws/2);
[N,Wn]=cheb1ord(Wp,Ws,Rp,Rs,'s');% 求原型模擬低通濾波器的階數和頻寬 
[bs,as]=cheby1(N,Rp,Wn,'s');     % 求模擬低通濾波器的係數
[B,A]=bilinear(bs,as,Fs);        % 按Fs把模擬低通濾波器的係數轉換成數字濾波器
[H1,f1]=freqz(B,A,1000,Fs);      % 計算數字濾波器的響應曲線
Wp=2/T*tan(wp/2);                % 把通帶和阻帶按Fs=1進行預畸
Ws=2/T*tan(ws/2);
[N,Wn]=cheb1ord(Wp,Ws,Rp,Rs,'s');% 求原型模擬低通濾波器的階數和頻寬 
[bs,as]=cheby1(N,Rp,Wn,'s');     % 求模擬低通濾波器的係數
[B,A]=bilinear(bs,as,1);         % 按Fs=1把模擬低通濾波器的係數轉換成數字濾波器
[H2,f2]=freqz(B,A,1000,Fs);      % 計算數字濾波器的響應曲線,恢復原取樣頻率,N為進行繪圖時的取樣點數,點數越大,濾波器頻率響應曲線越平滑
% 作圖
subplot 211; plot(f1,20*log10(abs(H1)),'k','linewidth',2)
xlabel('頻率/Hz'); ylabel('幅值/dB')
title('切比雪夫I型低通濾波器幅頻響應(bilinear中Fs=1000)')
axis([0 300 -50 5]); %grid; 
line([100 100],[-50 5],'color','k','linestyle',':');
line([200 200],[-50 5],'color','k','linestyle',':');
line([0 300],[-40 -40],'color','k','linestyle','--');
line([0 300],[-2 -2],'color','k','linestyle','--');
[H2,f2]=freqz(B,A,1000,Fs);
subplot 212; plot(f2,20*log10(abs(H2)),'k','linewidth',2)
xlabel('頻率/Hz'); ylabel('幅值/dB')
title('切比雪夫I型低通濾波器幅頻響應(bilinear中Fs=1)')
axis([0 300 -50 5]); %grid; 
line([100 100],[-50 5],'color','k','linestyle',':');
line([200 200],[-50 5],'color','k','linestyle',':');
line([0 300],[-40 -40],'color','k','linestyle','--');
line([0 300],[-2 -2],'color','k','linestyle','--');
set(gcf,'color','w') 

Fs取1或1000,效果是一樣的,因為在進行頻率畸變時候,對應的T=1/Fs也會不同。

4.直接設計數字濾波器

(1)函式呼叫

[bb,ab]=butter(N,Wn);        % 求巴特沃斯濾波器係數

不加‘s’就可以直接求出數字濾波器的係數

(2)例子(給出的實際頻率)

現設計一帶通濾波器,要求通帶頻率為1.5~10Hz,阻頻率為1和12Hz,Ap=3,As=15,設計巴特沃斯濾波器

clear all; clc; close all;
load bzsdata.mat                % 讀入資料
N=length(bzs);                  % 原始資料長
t=(0:N-1)/Fs;% 設定時間
fp1=[1.5 10]
fs1=[1 12]
wp1=2*fp1/Fs;
ws1=2*fs1/Fs;
Ap=3;
As=15;
[n,Wn]=buttord(wp1,ws1,Ap,As)
[b1,a1]=butter(n,Wn);
freqz(b1,a1,300,Fs)

然而設計出的濾波器頻率響應曲線為:
在這裡插入圖片描述
可以看到濾波器輸出產生了溢位和沒有數值。這是因為低通濾波器的截止頻率或帶通濾波器的中心頻率相對取樣頻率很低,或者濾波器的過渡帶很窄,則在設計濾波器後會發現濾波器階數比較大,這種情況下通常採用降低訊號的取樣頻率
在這裡插入圖片描述
現修改如下:

clear all; clc; close all;
load bzsdata.mat                % 讀入資料
N=length(bzs);                  % 原始資料長
t=(0:N-1)/Fs;                   % 設定時間
x=resample(bzs,1,5);            % 降取樣
N1=length(x);                   % 降取樣後的長度
fs=Fs/5;                        % 降取樣後的取樣頻率
fs2=fs/2;                       % 降取樣後取樣頻率的一半
t1=(0:N1-1)/fs;                 % 降取樣後的時間刻度
fp1=[1.5 10];                   % 通帶頻率
fs1=[1 12];                     % 阻帶頻率
wp1=fp1/fs2;                    % 歸一化通帶頻率
ws1=fs1/fs2;                    % 歸一化阻帶頻率
Ap=3; As=15;                    % 通帶波紋和阻帶衰減
[n,Wn]=buttord(wp1,ws1,Ap,As);  % 求濾波器原型階數和頻寬
[bn1,an1]=butter(n,Wn);         % 求數字濾波器係數
[H,f]=freqz(bn1,an1,1000,fs);   % 求數字濾波器幅頻曲線

y1=filter(bn1,an1,x);           % 對降取樣後的資料進行濾波
y=resample(y1,5,1);             % 對濾波器輸出恢復原取樣頻率
% 作圖
figure(1)
subplot 311; plot(t,bzs,'k');
xlabel('時間/秒'); title('原始資料波形')
subplot 312; plot(t1,x,'k');
xlabel('時間/秒'); title('降取樣後資料波形')
subplot 313; plot(t,y,'k')
xlabel('時間/秒'); title('濾波後資料波形')
set(gcf,'color','w');
figure(2)
plot(f,abs(H),'k');
grid; axis([0 25 0 1.1]); 
xlabel('頻率/Hz'); ylabel('幅值')
title('巴特沃斯濾波器的幅值響應')
set(gcf,'color','w');

注意,在實際應用中,buttord函式直接進行數字濾波器階數設計時wp,ws的計算方法如上段所示,直接進行歸一化,原頻率除以Fs/2。注意與模擬濾波器的區別

5.陷波器

(1)定義:

在訊號測量中,有時訊號會被一些干擾的正弦訊號所淹沒,所以要通過一些窄帶濾波器把這些正弦訊號濾除。窄帶濾波器的目的就是構成窄帶濾波,

6.全通濾波器

全通濾波器的應用之一是進行相位補償或相位均衡。具體說,由於IIR濾波器很難實現線性相位,因此在濾波器後加一個全通濾波器,調節全通濾波器的相位。

三、濾波使用函式

1.數字濾波

在這裡插入圖片描述

2.數字陷波器

(1)函式呼叫

[b,a]=iirnotch(Wo,Bw)

其中Wo為陷波器的中心頻率,Bw是陷波器的頻寬,引數b和a是數字濾波器的係數

3.全相位數字濾波器

(1)函式呼叫

[b,a]=iirgrpdelay(N,F,Edges,Gd)

其中N是全通濾波器階數,N必須是偶數,F和Gd是指在某一頻率區間內群延遲的指標,F和Gd都是向量,Gd是通帶的邊緣值,b和a是數字濾波器的係數。

四、濾波器的濾波過程用差分方程表示

差分方程實現濾波

相關文章