www.gusucode.com > 生物信号工具箱 > 生物信号工具箱/生物信号工具箱/Biosignal/Biosignal/Filters/reson.m
function [b,a]=reson(fn,fs,K,mthd) % % [b,a]=reson(fn,fs,K,mthd) % % Generates a resonant (single pole) filter.Parameters are: % % fn - desired location of pole (Hz) % fs - sampling frequency (Hz) % K - Q factor or pole magnitude(depending on witch method used % in the mthd parameter described below) % mthd - (string) desired method to use. Options are: % % 'poleMag' : in this case K specifies the magnitude of the pole in the filter % (0 < K <1). Notice that the gain might not be equal to 1 at the resonant % peak. And this method might provide *less* attenuaton at the sidebands than % method using 'Q2' below. % % 'Q1' : simple method utilizing a second order feedback system. K specificies % the filter's Q (how narrow the filter is). This method has significant % drawbacks, use Q2 instead. % % 'Q2' : improved version of Q1 with more symetrical response and unit gain at desired % pole location (see example). % % b - filter FIR taps % a - filter AR taps % % % %Example % fs=5000; % fn=50; % K=50; % r=0.997; % [b1,a1]=reson(fn,fs,K,'Q1'); % [b2,a2]=reson(fn,fs,K,'Q2'); % [b3,a3]=reson(fn,fs,r,'poleMag'); % % [H1,F]=freqz(b1,a1,500,fs); % hold on; % [H2,F]=freqz(b2,a2,500,fs); % [H3,F]= freqz(b3,a3,500,fs); % plot(F,20*log10(abs(H1)));hold on;grid on % plot(F,20*log10(abs(H2)),'r') % plot(F,20*log10(abs(H3)),'g') % legend('Q1','Q2','PoleMag') %Written By Ikaro Silva 2008 % % _______________________________________________________________________________ % Copyleft (l) 2010 by Ikaro Silva, All Rights Reserved % Contact Ikaro Silva (ikarosilva@ieee.org) % % This library (Biosignal Toolbox) is free software; you can redistribute % it and/or modify it under the terms of the GNU General Public License as published by % the Free Software Foundation; either version 2 of the License, or % (at your option) any later version. % % This program is distributed in the hope that it will be useful, % but WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details. % % You should have received a copy of the GNU General Public License % along with this program; if not, write to the Free Software % Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA % 02111-1307 USA % % _______________________________________________________________________________ switch mthd case 'poleMag' cW = cos(2*pi*fn/fs); b=[1]; a=[1 -2*K*cW K^2]; case 'Q1' % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Puts a pole exatly on the unit circle % at the fn Hz line. %Not recommended!!! Use Q2 instead cW =2*cos(2*pi*fn/fs); Hn=[1 -cW 1]; b=[1]; a=[1 K.*Hn]; case 'Q2' %Similar to 'Q1' but more accurate %For details see IEEE SP 2008 (5), pg 113 beta=1+K; f=pi*fn/fs; numA=tan(pi/4 - f); denA=sin(2*f)+cos(2*f)*numA; A=numA/denA; b=[1 -2*A A.^2]; a=[ (beta + K*(A^2)) -2*A*(beta+K) ((A^2)*beta + K)]; end