Matlab实现FIR滤波器的设计与应用

这个网址中详细介绍了关于matlab中fdatool模块的功能使用:https://wenku.baidu.com/view/3e6792de2af90242a995e558.html

一、FIR数字滤波器参数介绍

这个博文写的很详细↓↓↓↓
https://blog.csdn.net/zhoufan900428/article/details/8969470

主要是求Rs、Ap

1、只要是“衰减”:皆为-值,即阻带衰减40,除了矩形窗不能选,其余三种都可以
2、通过“过渡带”的公式求解N,即阶数。过渡带=ws-wp(=通带起始-阻带边缘)


二、FIR低通、高通、带通、带阻滤波器对应的代码

https://wenku.baidu.com/view/6c53300e79563c1ec5da71e8.html
里面用了汉宁窗、汉明窗、凯德窗。不过不是分别求,而是某一个特性滤波器只用了这一种方法。


三、利用fdatool生成滤波器

这是一个很取巧的方法,你只需要在里面输入对应的通带频率、阻带频率,它会自动生成一个滤波器。然后可以通过操作转换成matlab代码,该代码是以函数的形式出现的,里面已经将该滤波器做出来了,你只需要在filter()里面放入它的名字即可。
1、输入原始信号
clear;
clc;
[y,fs]=audioread('E:\son\dudu.wav'); %将声音放于matlab中
info=audioinfo('E:\son\dudu.wav')
%sound(y,fs);
T=1/fs; %采样时间
t=(0:length(y)-1)*T;%时间
f=(0:length(y)-1)*fs/length(y);
figure(1);
yz=y(:,1);%左声道
subplot(2,1,1);
plot(t,yz);%输入信号时域曲线
title('原始信号时域');
xlabel('时间');
ylabel('振幅');
subplot(2,1,2);
n=length(yz);%进行变换的点数
y1=fft(yz,n); %对n点进行傅里叶变换到频域
F=fs/length(yz); %谱分辨率,频谱间隔
plot(f,abs(y1));%左声道频谱图
axis([0,1000,0,10000]);
title('原始信号频谱');
xlabel('F(Hz)');
ylabel('H(jw)');
grid on
2、加噪声
noise1=sin(2*pi*100*t);
noise2=sin(2*pi*500*t);
noise3=sin(2*pi*900*t);
x1=noise1+noise2+noise3;
x1=yz+x1’;  #对x转置
figure(2);
subplot(2,1,1);
plot(t,x1);
title('加噪信号时域');
xlabel('时间');
ylabel('振幅');
subplot(2,1,2);
nx=length(x1);
x=fft(x1,nx);
plot(f,abs(x));%左声道频谱图
axis([0,1000,0,10000]);
title('加噪信号频谱');
xlabel('F(Hz)');
ylabel('H(jw)');
3、带阻滤波器100Hz
x2=filter(fir100,x1);
figure(3);
subplot(2,1,1);
plot(t,x2);
title('100HZ被滤掉后的时域');
xlabel('时间');
ylabel('振幅');
subplot(2,1,2);
plot(f,abs(fft(x2)));
axis([0 1000 0 10000]);
title('100Hz被滤后的频谱');
xlabel('F(Hz)');
ylabel('H(jw)');


function Hd = fir100
Fs = 48000;  #Sampling Frequency
Fpass1 = 90;             # First Passband Frequency
Fstop1 = 95;             # First Stopband Frequency
Fstop2 = 110;           # Second Stopband Frequency
Fpass2 = 120;           #Second Passband Frequency
Dpass1 = 0.028774368332;  # First Passband Ripple
Dstop  = 0.001;           # Stopband Attenuation
Dpass2 = 0.057501127785; #Second Passband Ripple
flag   = 'scale';         #Sampling Flag

#Calculate the order from the parameters using KAISERORD.
[N,Wn,BETA,TYPE] = kaiserord([Fpass1 Fstop1 Fstop2 Fpass2]/(Fs/2), [1 ...
                             0 1], [Dpass1 Dstop Dpass2]);

b  = fir1(N, Wn, TYPE, kaiser(N+1, BETA), flag);
Hd = dfilt.dffir(b);


4、带阻滤波器500Hz
x3=filter(fir500,x2);
figure(4);
subplot(2,1,1);
plot(t,x3);
title('100HZ、500HZ被滤掉后的时域');
xlabel('时间');
ylabel('振幅');
subplot(2,1,2);
plot(f,abs(fft(x3)));
axis([0 1000 0 10000]);
title('100Hz、500Hz被滤后的频谱');
xlabel('F(Hz)');
ylabel('H(jw)');


function Hd = fir500
 Fs = 48000;  # Sampling Frequency
 Fpass1 = 496;             3First Passband Frequency
 Fstop1 = 498;             #First Stopband Frequency
 Fstop2 = 504;           # Second Stopband Frequency
 Fpass2 = 508;             # Second Passband Frequency
 Dpass1 = 0.028774368332;  % First Passband Ripple
 Dstop  = 0.001;          # Stopband Attenuation
 Dpass2 = 0.057501127785;  % Second Passband Ripple
 flag   = 'scale';        # Sampling Flag

# Calculate the order from the parameters using KAISERORD.
 [N,Wn,BETA,TYPE] = kaiserord([Fpass1 Fstop1 Fstop2 
 Fpass2]/(Fs/2), [1 ...
                             0 1], [Dpass1 Dstop Dpass2]);

# Calculate the coefficients using the FIR1 function.
 b  = fir1(N, Wn, TYPE, kaiser(N+1, BETA), flag);
 Hd = dfilt.dffir(b);
5、带阻滤波器900Hz
x4=filter(fir900,x3);
figure(5);
subplot(2,1,1);
plot(t,x4);
title('全被滤掉后的时域');
xlabel('时间');
ylabel('振幅');
subplot(2,1,2);
plot(f,abs(fft(x4)));
axis([0 1000 0 10000]);
title('全部被滤后的频谱');
xlabel('F(Hz)');
ylabel('H(jw)');
sound(x4,fs);

function Hd = fir900
 Fpass1 = 890;             #First Passband Frequency
 Fstop1 = 895;             # First Stopband Frequency
 Fstop2 = 905;             # Second Stopband Frequency
 Fpass2 = 910;             # Second Passband Frequency
 Dpass1 = 0.028774368332;  # First Passband Ripple
 Dstop  = 0.001;           # Stopband Attenuation
 Dpass2 = 0.057501127785;  # Second Passband Ripple
 flag   = 'scale';         # Sampling Flag

# Calculate the order from the parameters using KAISERORD.
 [N,Wn,BETA,TYPE] = kaiserord([Fpass1 Fstop1 Fstop2 
 Fpass2]/(Fs/2), [1 ...
                             0 1], [Dpass1 Dstop Dpass2]);

#Calculate the coefficients using the FIR1 function.
 b  = fir1(N, Wn, TYPE, kaiser(N+1, BETA), flag);
 Hd = dfilt.dffir(b);
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 217,734评论 6 505
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 92,931评论 3 394
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 164,133评论 0 354
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 58,532评论 1 293
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 67,585评论 6 392
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 51,462评论 1 302
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 40,262评论 3 418
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 39,153评论 0 276
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 45,587评论 1 314
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,792评论 3 336
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,919评论 1 348
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,635评论 5 345
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 41,237评论 3 329
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,855评论 0 22
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,983评论 1 269
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 48,048评论 3 370
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,864评论 2 354

推荐阅读更多精彩内容