matlab下一个M文件里有用各个窗函数实现的滤波器,怎样调用这个滤波器

matlab下一个M文件里有用各个窗函数实现的滤波器,怎样调用这个滤波器,第1张

%本函数利用窗函数法设计带通滤波器,主要用来滤出单一频率,即中心频率

%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为文件,然后运行之。

欢迎分享,转载请注明来源:内存溢出

原文地址: http://outofmemory.cn/langs/12176705.html

(0)
打赏 微信扫一扫 微信扫一扫 支付宝扫一扫 支付宝扫一扫
上一篇 2023-05-21
下一篇 2023-05-21

发表评论

登录后才能评论

评论列表(0条)

保存