这个网址中详细介绍了关于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);