MATLAB辅助设计常用滤波器及其离散化
发布日期:2021-07-01 04:03:07
浏览次数:2
分类:技术文章
本文共 7080 字,大约阅读时间需要 23 分钟。
文章目录
声明:本文内容整理自网络,仅作为博主学习笔记记录,各部分版权归原作者所有。
几个概念理解:fdatool
工具箱,参考- 本文使用的方法
1. 三种滤波器的对比
传送门:
2. 巴特沃斯滤波器
测试用原始信号设置
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ButterWorthFilter.m% 巴特沃夫滤波器的设计%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear;close all;clc; fs = 1000; %Hz 采样频率Ts = 1/fs;N = 1000; %序列长度t = (0:N-1)*Ts;delta_f = 1*fs/N;f1 = 50;f2 = 100;f3 = 200;f4 = 400;x1 = 2*0.5*sin(2*pi*f1*t);x2 = 2*0.2*sin(2*pi*f2*t);x3 = 2*0.3*sin(2*pi*f3*t);x4 = 2*0.6*sin(2*pi*f4*t);x = x1 + x2 + x3 + x4; %待处理信号由四个分量组成 X = fftshift(abs(fft(x)))/N;X_angle = fftshift(angle(fft(x)));f = (-N/2:N/2-1)*delta_f; figure(1);subplot(3,1,1);plot(t,x);title('原信号');subplot(3,1,2);plot(f,X);grid on;title('原信号频谱幅度特性');subplot(3,1,3);plot(f,X_angle);title('原信号频谱相位特性');grid on;
(1)低通滤波器
[B,A] = butter(N,Wn,'low')
设计案例:
%设计一个巴特沃夫低通滤波器,要求把50Hz的频率分量保留,其他分量滤掉wp = 65/(fs/2); %通带截止频率,取50~100中间的值,并对其归一化ws = 85/(fs/2); %阻带截止频率,取50~100中间的值,并对其归一化alpha_p = 3; %通带允许最大衰减为 dbalpha_s = 20;%阻带允许最小衰减为 db%获取阶数和截止频率[ N1 wc1 ] = buttord( wp , ws , alpha_p , alpha_s);%获得转移函数系数[ b a ] = butter(N1,wc1,'low');%滤波filter_lp_s = filter(b,a,x);X_lp_s = fftshift(abs(fft(filter_lp_s)))/N;X_lp_s_angle = fftshift(angle(fft(filter_lp_s)));figure(2);freqz(b,a); %滤波器频谱特性figure(3);subplot(3,1,1);plot(t,filter_lp_s);grid on;title('低通滤波后时域图形');subplot(3,1,2);plot(f,X_lp_s);title('低通滤波后频域幅度特性');subplot(3,1,3);plot(f,X_lp_s_angle);title('低通滤波后频域相位特性');
(2)高通滤波器
[B,A] = butter(N,Wn,'high')
设计案例:
%设计一个高通滤波器,要求把400Hz的频率分量保留,其他分量滤掉wp = 350/(fs/2); %通带截止频率,取200~400中间的值,并对其归一化ws = 380/(fs/2); %阻带截止频率,取200~400中间的值,并对其归一化alpha_p = 3; %通带允许最大衰减为 dbalpha_s = 20;%阻带允许最小衰减为 db%获取阶数和截止频率[ N2 wc2 ] = buttord( wp , ws , alpha_p , alpha_s);%获得转移函数系数[ b a ] = butter(N2,wc2,'high');%滤波filter_hp_s = filter(b,a,x);X_hp_s = fftshift(abs(fft(filter_hp_s)))/N;X_hp_s_angle = fftshift(angle(fft(filter_hp_s)));figure(4);freqz(b,a); %滤波器频谱特性figure(5);subplot(3,1,1);plot(t,filter_hp_s);grid on;title('高通滤波后时域图形');subplot(3,1,2);plot(f,X_hp_s);title('高通滤波后频域幅度特性');subplot(3,1,3);plot(f,X_hp_s_angle);title('高通滤波后频域相位特性');
(3)带通滤波器
% [B,A] = butter(N,Wn)fs = 1000 ; %信号的采样频率wp=[8 30]*2/fs; %通带边界频率 ,单位为rad/sws=[7 32]*2/fs; %阻带边界频率 ,单位为rad/sRp=1; %通带最大波纹度 ,单位dB (不要太小)Rs=30; %表示阻带最小衰减,单位dB[N,Wn]=buttord(Wp,Ws,Rp,Rs); %巴特沃斯数字滤波器的阶数n和-3dB归一化截止频率Wn [B,A]=butter(N,Wn); %得到n阶巴特沃斯滤波的分子分母dataOut = filter(B,A,dataIn);
设计案例:
%设计一个带通滤波器,要求把50Hz和400Hz的频率分量滤掉,其他分量保留wp = [65 385 ] / (fs/2); %通带截止频率,50~100、200~400中间各取一个值,并对其归一化ws = [75 375 ] / (fs/2); %阻带截止频率,50~100、200~400中间各取一个值,并对其归一化alpha_p = 3; %通带允许最大衰减为 dbalpha_s = 20;%阻带允许最小衰减为 db%获取阶数和截止频率[ N3 wn ] = buttord( wp , ws , alpha_p , alpha_s);%获得转移函数系数[ b a ] = butter(N3,wn,'bandpass');%滤波filter_bp_s = filter(b,a,x);X_bp_s = fftshift(abs(fft(filter_bp_s)))/N;X_bp_s_angle = fftshift(angle(fft(filter_bp_s)));figure(6);freqz(b,a); %滤波器频谱特性figure(7);subplot(3,1,1);plot(t,filter_bp_s);grid on;title('带通滤波后时域图形');subplot(3,1,2);plot(f,X_bp_s);title('带通滤波后频域幅度特性');subplot(3,1,3);plot(f,X_bp_s_angle);title('带通滤波后频域相位特性');
(4)带阻滤波器
[B,A] = butter(N,Wn,'stop')
设计案例:
%设计一个带阻滤波器,要求把50Hz和400Hz的频率分量保留,其他分量滤掉wp = [65 385 ] / (fs/2); %通带截止频率?,50~100、200~400中间各取一个值,并对其归一化ws = [75 375 ] / (fs/2); %阻带截止频率?,50~100、200~400中间各取一个值,并对其归一化alpha_p = 3; %通带允许最大衰减为 dbalpha_s = 20;%阻带允许最小衰减为 db%获取阶数和截止频率[ N4 wn ] = buttord( wp , ws , alpha_p , alpha_s);%获得转移函数系数[ b a ] = butter(N4,wn,'stop');%滤波filter_bs_s = filter(b,a,x);X_bs_s = fftshift(abs(fft(filter_bs_s)))/N;X_bs_s_angle = fftshift(angle(fft(filter_bs_s)));figure(8);freqz(b,a); %滤波器频谱特性figure(9);subplot(3,1,1);plot(t,filter_bs_s);grid on;title('带阻滤波后时域图形');subplot(3,1,2);plot(f,X_bs_s);title('带阻滤波后频域幅度特性');subplot(3,1,3);plot(f,X_bs_s_angle);title('带阻滤波后频域相位特性');
3. 契比雪夫I型滤波器
(1)低通滤波器
fs = 1000; % Hz采样频率wp = 55/(fs/2); %通带截止频率,取50~100中间的值,并对其归一化ws = 90/(fs/2); %阻带截止频率,取50~100中间的值,并对其归一化Rp = 3; %通带允许最大衰减为 dbRs = 40; %阻带允许最小衰减为 db[n,Wn]=cheb1ord(Wp,Ws,Rp,Rs); % 获取阶数和截止频率[b,a]=cheby1(n,Rp,Wn, 'low'); %获得转移函数系数dataOut = filter(b,a,dataIn); %信号滤波运算
function y=lowp(x,f1,f3,rp,rs,Fs) %低通滤波 %使用注意事项:通带或阻带的截止频率的选取范围是不能超过采样率的一半 %即,f1,f3的值都要小于 Fs/2 %x:需要带通滤波的序列 % f 1:通带截止频率 % f 3:阻带截止频率 %rp:边带区衰减DB数设置 %rs:截止区衰减DB数设置 %FS:序列x的采样频率 % rp=0.1;rs=30;%通带边衰减DB值和阻带边衰减DB值 % Fs=2000;%采样率 % wp=2*pi*f1/Fs; ws=2*pi*f3/Fs; % 设计切比雪夫滤波器; [n,wn]=cheb1ord(wp/pi,ws/pi,rp,rs); [bz1,az1]=cheby1(n,rp,wp/pi); %查看设计滤波器的曲线 [h,w]=freqz(bz1,az1,256,Fs); h=20*log10(abs(h)); figure;plot(w,h);title('所设计滤波器的通带曲线');grid on; % y=filter(bz1,az1,x);%对序列x滤波后得到的序列yend
(2)高通滤波器
function y=highp(x,f1,f3,rp,rs,Fs) %高通滤波 %使用注意事项:通带或阻带的截止频率的选取范围是不能超过采样率的一半 %即,f1,f3的值都要小于 Fs/2 %x:需要带通滤波的序列 % f 1:通带截止频率 % f 2:阻带截止频率 %rp:边带区衰减DB数设置 %rs:截止区衰减DB数设置 %FS:序列x的采样频率 % rp=0.1;rs=30;%通带边衰减DB值和阻带边衰减DB值 % Fs=2000;%采样率 % wp=2*pi*f1/Fs; ws=2*pi*f3/Fs; % 设计切比雪夫滤波器; [n,wn]=cheb1ord(wp/pi,ws/pi,rp,rs); [bz1,az1]=cheby1(n,rp,wp/pi,'high'); %查看设计滤波器的曲线 [h,w]=freqz(bz1,az1,256,Fs); h=20*log10(abs(h)); figure;plot(w,h);title('所设计滤波器的通带曲线');grid on; y=filter(bz1,az1,x);end
(3)带通滤波器
function y=bandp(x,f1,f3,fsl,fsh,rp,rs,Fs) %带通滤波 %使用注意事项:通带或阻带的截止频率与采样率的选取范围是不能超过采样率的一半 %即,f1,f3,fs1,fsh,的值小于 Fs/2 %x:需要带通滤波的序列 % f 1:通带左边界 % f 3:通带右边界 % fs1:衰减截止左边界 % fsh:衰变截止右边界 %rp:边带区衰减DB数设置 %rs:截止区衰减DB数设置 %FS:序列x的采样频率 % f1=300;f3=500;%通带截止频率上下限 % fsl=200;fsh=600;%阻带截止频率上下限 % rp=0.1;rs=30;%通带边衰减DB值和阻带边衰减DB值 % Fs=2000;%采样率 % wp1=2*pi*f1/Fs; wp3=2*pi*f3/Fs; wsl=2*pi*fsl/Fs; wsh=2*pi*fsh/Fs; wp=[wp1 wp3]; ws=[wsl wsh]; % % 设计切比雪夫滤波器; [n,wn]=cheb1ord(ws/pi,wp/pi,rp,rs); [bz1,az1]=cheby1(n,rp,wp/pi); %查看设计滤波器的曲线 [h,w]=freqz(bz1,az1,256,Fs); h=20*log10(abs(h)); figure;plot(w,h);title('所设计滤波器的通带曲线');grid on; y=filter(bz1,az1,x);end
(4)带阻滤波器
function y=bands(x,f1,f3,fsl,fsh,rp,rs,Fs) %带阻滤波 %使用注意事项:通带或阻带的截止频率与采样率的选取范围是不能超过采样率的一半 %即,f1,f3,fs1,fsh,的值小于 Fs/2 %x:需要带通滤波的序列 % f 1:通带左边界 % f 3:通带右边界 % fs1:衰减截止左边界 % fsh:衰变截止右边界 %rp:边带区衰减DB数设置 %rs:截止区衰减DB数设置 %FS:序列x的采样频率 % f1=300;f3=500;%通带截止频率上下限 % fsl=200;fsh=600;%阻带截止频率上下限 % rp=0.1;rs=30;%通带边衰减DB值和阻带边衰减DB值 % Fs=2000;%采样率 % wp1=2*pi*f1/Fs; wp3=2*pi*f3/Fs; wsl=2*pi*fsl/Fs; wsh=2*pi*fsh/Fs; wp=[wp1 wp3]; ws=[wsl wsh]; % % 设计切比雪夫滤波器; [n,wn]=cheb1ord(ws/pi,wp/pi,rp,rs); [bz1,az1]=cheby1(n,rp,wp/pi,'stop'); %查看设计滤波器的曲线 [h,w]=freqz(bz1,az1,256,Fs); h=20*log10(abs(h)); figure;plot(w,h);title('所设计滤波器的通带曲线');grid on; y=filter(bz1,az1,x);end
4. 贝塞尔滤波器
% [b,a] = besself(n,Wn) %低通滤波器(Wn为二元向量时,为带通) % Wn是一个实际频率值,不是相对量[b,a] = besself(n,Wn,’ftype’); % ftype类型low,hight,stopdataOut = filter(b,a,dataIn); %信号滤波运算
5. MATLAB中Filter type类型表示
'low'
'bandpass'
'high'
'stop'
6. 滤波器的离散化
传送门:
参考文献:
- https://blog.csdn.net/xiaokun19870825/article/details/78777325
- https://blog.csdn.net/cqfdcw/article/details/84939698
- https://www.mathworks.com/help/signal/ref/cheby1.html
- https://www.mathworks.com/help/signal/ref/cheby1.html#bucqk89-ftype
- https://www.mathworks.com/help/matlab/ref/filter.html
- http://blog.sina.com.cn/s/blog_49c02a8c0100yszh.html
- https://blog.csdn.net/zhoufan900428/article/details/8969470
转载地址:https://miracle.blog.csdn.net/article/details/111212315 如侵犯您的版权,请留言回复原文章的地址,我们会给您删除此文章,给您带来不便请您谅解!
发表评论
最新留言
做的很好,不错不错
[***.243.131.199]2024年04月29日 00时55分17秒
关于作者
喝酒易醉,品茶养心,人生如梦,品茶悟道,何以解忧?唯有杜康!
-- 愿君每日到此一游!
推荐文章
node.js 实现一个简单的登录拦截器
2019-05-02
c++抽象类、纯虚函数以及巧用纯虚析构函数实现接口类【转】
2019-05-02
Caffe 安装错误记录及解决办法【转】
2019-05-02
Android用类继承Application的全局变量使用注意
2019-05-02
算法排序之桶排序
2019-05-02
lambda表达式初探
2019-05-02
C++ Template类模板的特化(3.3节, 3.4节)
2019-05-02
第05章 函数
2019-05-02
第08章 输入和输出
2019-05-02
QT中文乱码的解
2019-05-02
网上Qt多线程同步的一种普遍误识
2019-05-02
libcurl smtp发送邮件附件大小限制问题
2019-05-02
Qt中用QuaZip来压缩和解压缩文件
2019-05-02
第13章 Windows内存体系结构
2019-05-02
windows 和 linux 下c/c++内存分布(整理)
2019-05-02
Qt解析XML文件(QDomDocument)
2019-05-02
Qt图形视图框架
2019-05-02
Qt5中表格处理大数据量
2019-05-02
LeakCanary源码分析
2019-05-02
单例模式(Singleton)
2019-05-02