%本函数利用窗函数法设计带通滤波器,主要用来滤出单一频率,即中心频率
%data是输入的数据, centerFre是带通的中心频率, offsetFre是频偏,最终带通为centerFre +- offsetFre/2
%,sampFre是采样率
function y = BPassFilter(data, centerFre, offsetFre, sampFre)
%设计I型带通滤波器
M = 0 ; %滤波器阶数(必须是偶数)
Ap = 082; %通带衰减
As = 45; %阻带衰减
Wp1 = 2pi(centerFre - offsetFre)/sampFre; %算出下边频
Wp2 = 2pi(centerFre + offsetFre)/sampFre; %算出上边频
% (1)矩形窗
N = ceil(36sampFre/offsetFre); %计算滤波器阶数,采用矩形窗,3dB截频在中心频率到上下边频的中点
M = N - 1;
M = mod(M,2) + M ; %使滤波器为I型(偶数)
%单位脉冲响应的下脚标
h = zeros(1,M+1); %单位冲击响应变量赋初值
for k = 1:(M+1);
if (( k -1 - 05M)==0)
h(k) = Wp2/pi - Wp1/pi;
else
h(k) = Wp2sin(Wp2(k - 1 - 05M))/(pi(Wp2(k -1 - 05M))) - Wp1sin(Wp1(k - 1 - 05M))/(pi(Wp1(k -1 - 05M)));
end
end
% (2) Hann Window
% N = ceil(124sampFre/offsetFre); %计算滤波器阶数,采用矩形窗,3dB截频在中心频率到上下边频的中点
% M = N - 1;
% M = mod(M,2) + M ; %使滤波器为I型(偶数)
% h = zeros(1,M+1); %单位冲击响应变量赋初值
% for k = 1:(M+1);
% if (( k -1 - 05M)==0)
% h(k) = Wp2/pi - Wp1/pi;
% else
% h(k) = Wp2sin(Wp2(k - 1 - 05M))/(pi(Wp2(k -1 - 05M))) - Wp1sin(Wp1(k - 1 - 05M))/(pi(Wp1(k -1 - 05M)));
% end
% end
% K = 0:M;
% w = 05 - 05cos(2piK/M);
% h = hw;
% (3)Hamming Window
% N = ceil(14sampFre/offsetFre); %计算滤波器阶数,采用矩形窗,3dB截频在中心频率到上下边频的中点
% M = N - 1;
% M = mod(M,2) + M ; %使滤波器为I型(偶数)
% h = zeros(1,M+1); %单位冲击响应变量赋初值
% for k = 1:(M+1);
% if (( k -1 - 05M)==0)
% h(k) = Wp2/pi - Wp1/pi;
% else
% h(k) = Wp2sin(Wp2(k - 1 - 05M))/(pi(Wp2(k -1 - 05M))) - Wp1sin(Wp1(k - 1 - 05M))/(pi(Wp1(k -1 - 05M)));
% end
% end
% K = 0:M;
% w = 054 - 046cos(2pik/M);
% h = hw;
% (4)Blackman window
% N = ceil(228sampFre/offsetFre); %计算滤波器阶数,采用矩形窗,3dB截频在中心频率到上下边频的中点
% M = N - 1;
% M = mod(M,2) + M ; %使滤波器为I型(偶数)
% h = zeros(1,M+1); %单位冲击响应变量赋初值
% for k = 1:(M+1);
% if (( k -1 - 05M)==0)
% h(k) = Wp2/pi - Wp1/pi;
% else
% h(k) = Wp2sin(Wp2(k - 1 - 05M))/(pi(Wp2(k -1 - 05M))) - Wp1sin(Wp1(k - 1 - 05M))/(pi(Wp1(k -1 - 05M)));
% end
% end
% K = 0:M;
% w = 042 - 05cos(2piK/M) + 008cos(4piK/M);
% h = hw;
y = filter(h,[1],data);
SciPy提供了firwin用窗函数设计低通滤波器,firwin的调用形式如下:
firwin(N, cutoff, width=None, window='hamming')
其中N为滤波器的长度;cutoff为以正规化的频率;window为所使用的窗函数。
首先理想滤波器是不可能实现的,只能尽可能逼近理想滤波器的特性,但是由于成本和技术的复杂性的限制,在实际设计中就允许通带和阻带中存在一定的误差,即通带不一定是完全水平的,阻带信号也不一定都衰减到零,在通带与阻带之间还有一定宽度的过渡带 这些要求决定了滤波器的复杂程度,一般来说要求波动越小,就越难实现,而且这些因素也是相互制约的,鱼和熊掌不可兼得, 比如通带平整就意味着过渡带宽,衰减缓慢, 要求过渡带窄通带波动就大等等,故需要综合考虑,合理选择
以选择窗口函数为例, 要熟悉各种窗口函数的特性,根据具体要求来选择
最小阻带衰减 过渡带带宽△w
矩形窗 209dB 092π/M
汉宁窗 439dB 311π/M
海明窗 545dB 332π/M
布莱克曼窗 753dB 556π/M
将模拟频率转化为数字频率,设取样时间为T(要满足抽样定理) Ωp=2πfpT Ωs=2πfsT 过渡带宽度△Ω=Ωp-Ωs 阻带衰减已经超过74db,要选用Kaiser窗了,Kaiser的参数可变,要根据公式确定滤波器的参数一般都选用Ⅰ型线性相位滤波器即滤波器阶数M为偶数,程序如下: wp=;ws=;Ap=1;As=100; Rp=1-10^(-005Ap);Rs=10^(-005As); f=[fp fs]; a=[0 1]; dev=[Rp Rs]; [M,wc,beta,ftype]=kaiserord(f,a,dev); M=mod(M,2)+M; h=fir1(M,wc,ftype,kaiser(M+1,beta)); omega=linspace(0,pi,512); mag=freqz(h,[1],omega); plot(omega/pi,20log10(abs(mag))); grid; omega1=linspace(0,wp,512); h1=freqz(h,[1],omega1); omega2=linspace(ws,pi,512); h2=freqz(h,[1],omega2); fprintf('Ap=%4f\n',-20log10(min(abs(h1)))); fprintf('As=%4f\n',-20log10(max(abs(h2)))); 运行程序可以得到滤波器的通阻带衰减,画出频率响应,若同阻带衰减不满足要求还可以使用滤波器的优化,一般使用的等波纹FIR进行优化
上面这对代码分两部分:
1 %理想低通滤波器单位冲激响应函数
function hd=ideal_lp1(wc,N) % 这一行去掉分号
pha=angle(H);
2 % 主程序
clear all;
axis([0,1,-100,10])
其中1部分保存成一个叫ideal_lp1m的文件,放好别动;2部分保存成任意名字的m为文件,然后运行之。
欢迎分享,转载请注明来源:内存溢出
评论列表(0条)