% IIR濾波器設計
% 目的：設計一個採樣頻率為1000Hz、通帶截止頻率為50Hz、阻帶截止頻率為100Hz的低通濾波器，並要求通帶最大衰減為1dB，阻帶最小衰減為60dB。
clc;clear;close all;
% 1. 產生信號(頻率為10Hz和100Hz的正弦波疊加)
Fs=1000; % 採樣頻率1000Hz
t=0:1/Fs:1;
s10=sin(20*pi*t); % 產生10Hz正弦波
s100=sin(200*pi*t); % 產生100Hz正弦波
s=s10+s100; % 信號疊加
figure(1); % 畫圖
subplot(2,1,1);plot(s);grid;
title('原始信號');
% 2. FFT分析信號頻譜
len = 512;
y=fft(s,len); % 對信號做len點FFT變換
f=Fs*(0:len/2 - 1)/len;
subplot(2,1,2);plot(f,abs(y(1:len/2)));grid;
title('原始信號頻譜')
xlabel('Hz');ylabel('幅值');
% 3. IIR濾波器設計
N=0; % 階數
Fp=50; % 通帶截止頻率50Hz
Fc=100; % 阻帶截止頻率100Hz
Rp=1; % 通帶波紋最大衰減為1dB
Rs=60; % 阻帶衰減為60dB
% 3.0 計算最小濾波器階數
na=sqrt(10^(0.1*Rp)-1);
ea=sqrt(10^(0.1*Rs)-1);
N=ceil(log10(ea/na)/log10(Fc/Fp));
% 3.1 巴特沃斯濾波器
Wn=Fp*2/Fs;
[Bb Ba]=butter(N,Wn,'low'); % 調用MATLAB butter函數快速設計濾波器
[BH,BW]=freqz(Bb,Ba); % 繪製頻率響應曲線
Bf=filter(Bb,Ba,s); % 進行低通濾波
By=fft(Bf,len); % 對信號f1做len點FFT變換
figure(2); % 畫圖
subplot(2,1,1);plot(t,s10,'blue',t,Bf,'red');grid;
legend('10Hz原始信號','巴特沃斯濾波器濾波後');
subplot(2,1,2);plot(f,abs(By(1:len/2)));grid;
title('巴特沃斯低通濾波後信號頻譜');
xlabel('Hz');ylabel('幅值');
% 3.2 切比雪夫I型濾波器
[C1b C1a]=cheby1(N,Rp,Wn,'low'); % 調用MATLAB cheby1函數快速設計低通濾波器
[C1H,C1W]=freqz(C1b,C1a); % 繪製頻率響應曲線
C1f=filter(C1b,C1a,s); % 進行低通濾波
C1y=fft(C1f,len); % 對濾波後信號做len點FFT變換
figure(3); % 畫圖
subplot(2,1,1);plot(t,s10,'blue',t,C1f,'red');grid;
legend('10Hz原始信號','切比雪夫I型濾波器濾波後');
subplot(2,1,2);plot(f,abs(C1y(1:len/2)));grid;
title('切比雪夫I型濾波後信號頻譜');
xlabel('Hz');ylabel('幅值');
% 3.3 切比雪夫II型濾波器
[C2b C2a]=cheby2(N,Rs,Wn,'low'); % 調用MATLAB cheby2函數快速設計低通濾波器
[C2H,C2W]=freqz(C2b,C2a); % 繪製頻率響應曲線
C2f=filter(C2b,C2a,s); % 進行低通濾波
C2y=fft(C2f,len); % 對濾波後信號做len點FFT變換
figure(4); % 畫圖
subplot(2,1,1);plot(t,s10,'blue',t,C2f,'red');grid;
legend('10Hz原始信號','切比雪夫II型濾波器濾波後');
subplot(2,1,2);plot(f,abs(C2y(1:len/2)));grid;
title('切比雪夫II型濾波後信號頻譜');
xlabel('Hz');ylabel('幅值');
% 3.4 橢圓濾波器
[Eb Ea]=ellip(N,Rp,Rs,Wn,'low'); % 調用MATLAB ellip函數快速設計低通濾波器
[EH,EW]=freqz(Eb,Ea); % 繪製頻率響應曲線
Ef=filter(Eb,Ea,s); % 進行低通濾波
Ey=fft(Ef,len); % 對濾波後信號做len點FFT變換
figure(5); % 畫圖
subplot(2,1,1);plot(t,s10,'blue',t,Ef,'red');grid;
legend('10Hz原始信號','橢圓濾波器濾波後');
subplot(2,1,2);plot(f,abs(Ey(1:len/2)));grid;
title('橢圓濾波後信號頻譜');
xlabel('Hz');ylabel('幅值');
% 3.5 yulewalk濾波器
fyule=[0 Wn Fc*2/Fs 1]; % 在此進行的是簡單設計，實際需要多次仿真取最佳值
myule=[1 1 0 0]; % 在此進行的是簡單設計，實際需要多次仿真取最佳值
[Yb Ya]=yulewalk(N,fyule,myule); % 調用MATLAB yulewalk函數快速設計低通濾波器
[YH,YW]=freqz(Yb,Ya); % 繪製頻率響應曲線
Yf=filter(Yb,Ya,s); % 進行低通濾波
Yy=fft(Yf,len); % 對濾波後信號做len點FFT變換
figure(6); % 畫圖
subplot(2,1,1);plot(t,s10,'blue',t,Yf,'red');grid;
legend('10Hz原始信號','yulewalk濾波器濾波後');
subplot(2,1,2);plot(f,abs(Yy(1:len/2)));grid;
title('yulewalk濾波後信號頻譜');
xlabel('Hz');ylabel('幅值');
% 4. 各個濾波器的幅頻響應對比分析
A1 = BW*Fs/(2*pi);
A2 = C1W*Fs/(2*pi);
A3 = C2W*Fs/(2*pi);
A4 = EW*Fs/(2*pi);
A5 = YW*Fs/(2*pi);
figure(7); % 畫圖
subplot(1,1,1);plot(A1,abs(BH),A2,abs(C1H),A3,abs(C2H),A4,abs(EH),A5,abs(YH));grid;
xlabel('頻率／Hz');
ylabel('頻率響應幅度');
legend('butter','cheby1','cheby2','ellip','yulewalk');
% IIR濾波器設計
% 目的：設計一個採樣頻率為1000Hz、通帶截止頻率為50Hz、阻帶截止頻率為100Hz的低通濾波器，並要求通帶最大衰減為1dB，阻帶最小衰減為60dB。
clc;clear;close all;
% 1. 產生信號(頻率為10Hz和100Hz的正弦波疊加)
Fs=1000; % 採樣頻率1000Hz
t=0:1/Fs:1;
s10=sin(20*pi*t); % 產生10Hz正弦波
s100=sin(200*pi*t); % 產生100Hz正弦波
s=s10+s100; % 信號疊加
figure(1); % 畫圖
subplot(2,1,1);plot(s);grid;
title('原始信號');
% 2. FFT分析信號頻譜
len = 512;
y=fft(s,len); % 對信號做len點FFT變換
f=Fs*(0:len/2 - 1)/len;
subplot(2,1,2);plot(f,abs(y(1:len/2)));grid;
title('原始信號頻譜')
xlabel('Hz');ylabel('幅值');
% 3. IIR濾波器設計
N=0; % 階數
Fp=50; % 通帶截止頻率50Hz
Fc=100; % 阻帶截止頻率100Hz
Rp=1; % 通帶波紋最大衰減為1dB
Rs=60; % 阻帶衰減為60dB
% 3.0 計算最小濾波器階數
na=sqrt(10^(0.1*Rp)-1);
ea=sqrt(10^(0.1*Rs)-1);
N=ceil(log10(ea/na)/log10(Fc/Fp));
% 3.1 巴特沃斯濾波器
Wn=Fp*2/Fs;
[Bb Ba]=butter(N,Wn,'low'); % 調用MATLAB butter函數快速設計濾波器
[BH,BW]=freqz(Bb,Ba); % 繪製頻率響應曲線
Bf=filter(Bb,Ba,s); % 進行低通濾波
By=fft(Bf,len); % 對信號f1做len點FFT變換
figure(2); % 畫圖
subplot(2,1,1);plot(t,s10,'blue',t,Bf,'red');grid;
legend('10Hz原始信號','巴特沃斯濾波器濾波後');
subplot(2,1,2);plot(f,abs(By(1:len/2)));grid;
title('巴特沃斯低通濾波後信號頻譜');
xlabel('Hz');ylabel('幅值');
% 3.2 切比雪夫I型濾波器
[C1b C1a]=cheby1(N,Rp,Wn,'low'); % 調用MATLAB cheby1函數快速設計低通濾波器
[C1H,C1W]=freqz(C1b,C1a); % 繪製頻率響應曲線
C1f=filter(C1b,C1a,s); % 進行低通濾波
C1y=fft(C1f,len); % 對濾波後信號做len點FFT變換
figure(3); % 畫圖
subplot(2,1,1);plot(t,s10,'blue',t,C1f,'red');grid;
legend('10Hz原始信號','切比雪夫I型濾波器濾波後');
subplot(2,1,2);plot(f,abs(C1y(1:len/2)));grid;
title('切比雪夫I型濾波後信號頻譜');
xlabel('Hz');ylabel('幅值');
% 3.3 切比雪夫II型濾波器
[C2b C2a]=cheby2(N,Rs,Wn,'low'); % 調用MATLAB cheby2函數快速設計低通濾波器
[C2H,C2W]=freqz(C2b,C2a); % 繪製頻率響應曲線
C2f=filter(C2b,C2a,s); % 進行低通濾波
C2y=fft(C2f,len); % 對濾波後信號做len點FFT變換
figure(4); % 畫圖
subplot(2,1,1);plot(t,s10,'blue',t,C2f,'red');grid;
legend('10Hz原始信號','切比雪夫II型濾波器濾波後');
subplot(2,1,2);plot(f,abs(C2y(1:len/2)));grid;
title('切比雪夫II型濾波後信號頻譜');
xlabel('Hz');ylabel('幅值');
% 3.4 橢圓濾波器
[Eb Ea]=ellip(N,Rp,Rs,Wn,'low'); % 調用MATLAB ellip函數快速設計低通濾波器
[EH,EW]=freqz(Eb,Ea); % 繪製頻率響應曲線
Ef=filter(Eb,Ea,s); % 進行低通濾波
Ey=fft(Ef,len); % 對濾波後信號做len點FFT變換
figure(5); % 畫圖
subplot(2,1,1);plot(t,s10,'blue',t,Ef,'red');grid;
legend('10Hz原始信號','橢圓濾波器濾波後');
subplot(2,1,2);plot(f,abs(Ey(1:len/2)));grid;
title('橢圓濾波後信號頻譜');
xlabel('Hz');ylabel('幅值');
% 3.5 yulewalk濾波器
fyule=[0 Wn Fc*2/Fs 1]; % 在此進行的是簡單設計，實際需要多次仿真取最佳值
myule=[1 1 0 0]; % 在此進行的是簡單設計，實際需要多次仿真取最佳值
[Yb Ya]=yulewalk(N,fyule,myule); % 調用MATLAB yulewalk函數快速設計低通濾波器
[YH,YW]=freqz(Yb,Ya); % 繪製頻率響應曲線
Yf=filter(Yb,Ya,s); % 進行低通濾波
Yy=fft(Yf,len); % 對濾波後信號做len點FFT變換
figure(6); % 畫圖
subplot(2,1,1);plot(t,s10,'blue',t,Yf,'red');grid;
legend('10Hz原始信號','yulewalk濾波器濾波後');
subplot(2,1,2);plot(f,abs(Yy(1:len/2)));grid;
title('yulewalk濾波後信號頻譜');
xlabel('Hz');ylabel('幅值');
% 4. 各個濾波器的幅頻響應對比分析
A1 = BW*Fs/(2*pi);
A2 = C1W*Fs/(2*pi);
A3 = C2W*Fs/(2*pi);
A4 = EW*Fs/(2*pi);
A5 = YW*Fs/(2*pi);
figure(7); % 畫圖
subplot(1,1,1);plot(A1,abs(BH),A2,abs(C1H),A3,abs(C2H),A4,abs(EH),A5,abs(YH));grid;
xlabel('頻率／Hz');
ylabel('頻率響應幅度');
legend('butter','cheby1','cheby2','ellip','yulewalk');
效果圖
轉載請註明文章來源 – http://blog.csdn.net/v_hyx ，請勿用於任何商業用途
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
