MATLAB訊號處理——數字濾波器的設計
一、模擬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是數字濾波器的係數。
四、濾波器的濾波過程用差分方程表示
相關文章
- 語音訊號預處理——數字濾波器音訊
- 數字訊號處理實驗(四):數字濾波器結構
- 數字訊號處理基礎----插值、抽取濾波器
- 數字訊號處理:運用FFT簡單濾波FFT
- FPGA數字訊號處理(24)數字相關器設計(簡化結構)FPGA
- FPGA數字訊號處理(25)數字相關器設計(經典結構)FPGA
- FPGA數字訊號處理(26)加擾器與解擾器設計FPGA
- 數字濾波器和模擬濾波器(一)
- 數字訊號處理實驗一(離散時間訊號的MATLAB實現)Matlab
- 演算法 | 數字影像處理之「中值濾波」演算法
- MATLAB數字訊號處理(2)LFM脈衝雷達回波處理模擬Matlab
- FPGA數字訊號處理(27)卷積編碼器與Viterbi譯碼器設計FPGA卷積Viterbi
- MATLAB數字訊號處理(1)四種經典功率譜估計方法比較Matlab
- FPGA數字訊號處理(23)FSK解調技術(包絡檢波法)FPGA
- MATLAB及其訊號處理基礎Matlab
- matlab 濾波器中用到的函式Matlab函式
- 大牛講解訊號與系統以及數字訊號處理
- OpenCV計算機視覺學習(4)——影像平滑處理(均值濾波,高斯濾波,中值濾波,雙邊濾波)OpenCV計算機視覺
- 數字訊號處理c語言程式集C語言
- [00]數字影像處理-matlab速成Matlab
- [Python影象處理] 四.影象平滑之均值濾波、方框濾波、高斯濾波及中值濾波Python
- 訊號處理基本引數
- 音訊降噪-fir濾波器音訊
- 數字影像處理-取樣量化(Matlab)Matlab
- 高通WCD9375音訊編解碼器/數字濾波器晶片簡介音訊晶片
- FPGA數字訊號處理(22)FSK調製技術FPGA
- 資料平滑處理-均值|中值|Savitzky-Golay濾波器Go
- Python 影像處理 OpenCV (7):影像平滑(濾波)處理PythonOpenCV
- MATLAB數字影象處理(二)影象增強Matlab
- C++數字訊號處理演算法庫SP++C++演算法
- 利用Matlab filterDesigner 工具生成FIR濾波器函式,並呼叫實現低通濾波MatlabFilter函式
- 數字訊號處理:線性卷積、迴圈卷積、圓周卷積計算卷積
- 影像處理技術(二)濾波去噪(上)
- FPGA數字訊號處理(十三)鎖相環位同步技術的實現FPGA
- matlab生成常用訊號(方波、三角波、隨機訊號、單位衝激)Matlab隨機
- 詳解數字影像的濾波和邊緣檢測
- 點雲濾波器與過濾器過濾器
- 小波變換檢測訊號突變點的MATLAB實現Matlab