有些技术含量用法。filter本身作用是求卷积和conv
filter(B,1,X,[],dim);dim缺省为1,是按列滤波的,如果改为2,则是按行滤波。 y = filter(b, a, x),其中b,a为滤波器系数。计算系统在输入x作用下的零状态响应 y[k] 举例:计算低通滤波器的冲激响应
例题1 点平均滤波
f1=3;f2=40;fs=100;t=0:1/fs:1;
x=sin(2*pi*t*f1)+.25*sin(2*pi*t*f2);
b=ones(1,10)/10; y=filter(b,1,x);求冲激响应 stem(y);
yy=filtfilt(b,1,x); plot(t,x);hold on; plot(t,x,'r',t,yy,'g')
例 2利用 filter函数求滑动平均
Matlab有多种计算滑动平均的方法,现介绍基于 filter函数的计算方法。设原始数据为x, 平均窗口设为a(a为正整数),那么无权重滑动平均后的数据 y为: windowSize = a;
y=filter(ones(1,windowSize)/windowSize,1,x); 上述命令实际上计算的是: y(1)=(1/a)*x(1);
y(2)=(1/a)*x(2)+(1/a)*x(1);
y(a)=(1/a)*x(a)+(1/a)*x(a-1)+...+(1/a)*x(1); y(i)=(1/a)*x(i)+(1/a)*x(i-1)+...+(1/a)*x(i-a+1);
4. frezq数字滤波器的频率响应
[H, W]=freqz(B, A, N) 当 N是一个整数时返回 N点的频率向量H和N点的幅频响 应向量W。N最好选用 2的整数次幂,这样便于使用 FFT进行快速算法。
H为滤波器的复数放大倍数,w为频率向量,如果只想获得放大倍数的幅值,可以用 plot(w,abs(h)).如果是滤波器设计 plot(w/pi,abs(h))
滤波器放大倍数,低频时为1,高频时为0,即低通滤波器
N个频率点均匀地分布在单位圆的上半圆上。如果 N没有确定则却缺省为 512个点。 freqz(B, A, N) 将直接绘制频率响应图,而不返回任何值。
[H, W]=freqz(B, A, N, 'whole') 运用分布在整个单位圆上的 N个点。
H=freqz(B, A, W) 返回指定在 W向量中频率范围的频率响应,其中 W是以弧度为单位 在[0, pi]范围内。
[H,F]=freqz(B, A, N, Fs), [H,F]=freqz(B, A, N, Fs) 这两个函数给出了采样 频率Fs,则返回频率向量F,它们的单位都是Hz。
invfreqz()是其逆函数,它运用最小二乘法从已知的频率响应中求出传递函数模型。
例子 2
滤波器的特性分析二
3.数字滤波器的冲击响应:impz()
[H, T]=impz(B, A) 返回滤波器的冲击响应列向量 H和时间即采样间隔列向量T。 滤波器是用传递函数模型来限定的。
impz(B, A) 将直接绘制滤波器的冲击响应图。
例:设计一个高通 Chebyshev2型滤波器,它的具体要求是阻带的截止频率为10Hz,通带 的截止频率为30Hz,在通带内的最大衰减不超过0.1dB,在阻带内的最小衰减不小于40dB。 程序设计如下:
wp=30; ws=10; rp=0.1; rs=40; Fs=100; wp=10*2*pi; ws=30*2*pi;
[n, wn]=cheb2ord(wp/(Fs/2), ws/(Fs/2), rp, rs, 's'); [z, p, k]=cheb2ap(n, rs);
wn=wn*(Fs/2)*2*pi;
[A, B, C, D]=zp2ss(z,p,k);
[AT, BT, CT, DT]=lp2hp(A, B, C, D, wn);
[AT1, BT1, CT1, DT1]=bilinear(AT, BT, CT, DT, Fs); [num, den]=ss2tf(AT1, BT1, CT1, DT1); [H, W]=freqz(num, den); figure;
plot(W*Fs/(2*pi), abs(H)); grid; xlabel('频率/Hz'); ylabel('幅值');
impz(num, den);
一、巴特沃斯 IIR滤波器的设计
1 首先根据滤波器设计要求用 buttord函数求出最小滤波器阶数和截止频率
[n,Wn]= buttord(Wp,Ws,Rp,Rs)
其中Wp和 Ws分别是通带和阻带的拐角频率(截止频率),其取值范围为0至1之间。当 其值为 1时代表采样频率的一半。Rp和 Rs分别是通带和阻带区的波纹系数。 不同类型(高通、低通、带通和带阻)滤波器对应的 Wp和 Ws值遵循以下规则:
1.高通滤波器:Wp和 Ws为一元矢量且Wp>Ws; 2.低通滤波器:Wp和 Ws为一元矢量且Wp 3.带通滤波器:Wp和 Ws为二元矢量且Wp Butter函数可设计低通、高通、带通和带阻的数字和模拟 IIR滤波器,其特性为使通带内 的幅度响应最大限度地平坦,但同时损失截止频率处的下降斜度。在期望通带平滑的情况 下,可使用 butter函数。 butter函数的用法为: [b,a]=butter(n,Wn,ftype)。其中 n代表滤波器阶数,Wn代表滤波器的截止频率 二、契比雪夫 I型 IIR滤波器的设计 在期望通带下降斜率大的场合,应使用椭圆滤波器或契比雪夫滤波器。在 MATLAB下可使用 cheby1函数设计出契比雪夫 I型 IIR滤波器。 cheby1函数可设计低通、高通、带通和带阻契比雪夫 I型滤 IIR波器,其通带内为等波纹, 阻带内为单调。契比雪夫 I型的下降斜度比 II型大,但其代价是通带内波纹较大。 cheby1函数的用法为: [b,a]=cheby1(n,Rp,Wn,/ftype/) 在使用 cheby1函数设计 IIR滤波器之前,可使用 cheblord函数求出滤波器阶数 n和截止 频率Wn。cheblord函数可在给定滤波器性能的情况下,选择契比雪夫 I型滤波器的最小阶 和截止频率Wn。 cheblord函数的用法为: [n,Wn]=cheblord(Wp,Ws,Rp,Rs) 其中 Wp和 Ws分别是通带和阻带的拐角频率(截止频率),其取值范围为0至1之间。当 其值为 1时代表采样频率的一半。Rp和 Rs分别是通带和阻带区的波纹系数。 例1 选择设计 IIR的 Butterworth低通滤波器,其 Fs=22050Hz,Fp1=3400Hz, Fs1=5000Hz,Rp=2dB,Rs=20dB Fs=22050;Fp1=3400;Fs1=5000;Rp=3;Rs=20;%设计指标 wp1=2*Fp1 /Fs;ws1=2*Fs1 /Fs;%求归一化频率 % 确定 butterworth 的最小阶数 N 和频率参数 Wn [n,Wn]=buttord(wp1,ws1,Rp,Rs); [B,A] = butter(N,Wn);%确定传递函数的分子、分母系数 [h,f]=freqz(b,a,Nn,Fs_value);%生成频率响应参数 plot(f,20*log(abs(h))) %画幅频响应图 plot(f,angle(h)); %画相频响应图 %[N, Wn] = buttord(Wp, Ws, Rp, Rs) 确定 butterworth 的 N 和 Wn %[N, Wn] = cheblord ( (Wp, Ws, Rp, Rs) 确定 Chebyshev滤波器的 N 和 Wn %[N, Wn] = cheb2ord (Wp, Ws, Rp, Rs) 确定 Chebyshev2滤波器的 N 和 Wn %[N, Wn] = ellipord (Wp, Ws, Rp, Rs)确定椭圆(Ellipse) 滤波器的 N 和 Wn %[B,A] = butter(N,Wn,'type')设计'type'型巴特沃斯(Butterworth)滤波器 filter. %[B,A] = cheby1 (N,R,Wn, 'type')设计'type'型切比雪夫Ⅰ滤波器 filter. %[B,A] = cheby2(N,R,Wn, 'type')设计'type'型切比雪夫Ⅱ滤波器 filter. %[B,A] = ellip(N,Rp,Rs,Wn, 'type')设计'type' 型椭圆 filter. 例子 2 %实现了对频率为 20和 200Hz单频叠加 cos信号的低通滤波,使输出仅含有 20Hz分量 clear; fs=1200; %采样频率 N=1200; % N/fs 秒数据 n=0:N-1; t=n/fs; %时间 fL=20; fH=200; sL=cos(2*pi*fL*t); sH=cos(2*pi*fH*t); s=sL+sH; % s_in_f=fft(s); % i=1:250; % plot(i,s_in_f(i)); figure(1); plot(t,s); title('输入信号'); xlabel('t/s'); ylabel('幅度'); %设计低通滤波器: Wp = 50/fs; Ws = 100/fs; %截止频率50Hz,阻带截止频率100Hz,采样频率 200Hz [n,Wn] = buttord(Wp,Ws,1,50); %阻带衰减大于50db,通带纹波小于 1db %估算得到 Butterworth低通滤波器的最小阶数 N和 3dB截止频率 Wn [a,b]=butter(n,Wn); %设计 Butterworth低通滤波器 [h,f]=freqz(a,b,'whole',fs); %求数字低通滤波器的频率响应 f=(0:length(f)-1)'*fs/length(f); %进行对应的频率转换 figure(2); plot(f,abs(h)); %绘制 Butterworth低通滤波器的幅频响应图 title('巴氏低通滤波器'); grid; sF=filter(a,b,s); %叠加函数 s经过低通滤波器以后的新函数 figure(3); plot(t,sF); %绘制叠加函数 S经过低通滤波器以后的时域图形 xlabel('t/s'); ylabel('幅度'); SF=fft(sF,N); %对叠加函数S经过低通滤波器以后的新函数进行256点的基—2快 速傅立叶变换 mag=abs(SF); %求幅值 f=(0:length(SF)-1)'*fs/length(SF); %进行对应的频率转换 figure(4); plot(f,mag); %绘制叠加函数 S经过低通滤波器以后的频谱图 title('低通滤波后的频谱图'); 4 窗函数法 FIR滤波器设计实验 FIR滤波器的设计任务是选择有限长度的h(n)。使传输函数H( )满足技术要求。FIR滤波 器的设计方法有多种,如窗函数法、频率采样法及其它各种优化设计方法,本实验介绍窗 函数法的 FIR滤波器设计。 窗函数法是使用矩形窗、三角窗、巴特利特窗、汉明窗、汉宁窗和布莱克曼窗等设计出标 准响应的高通、低通、带通和带阻 FIR滤波器。 一、firl函数的使用 在 MATLAB下设计标准响应 FIR滤波器可使用 firl函数。firl函数以经典方法实现加窗线 性相位 FIR滤波器设计,它可以设计出标准的低通、带通、高通和带阻滤波器。firl函数 的用法为: b=firl(n,Wn,/ftype/,Window) 各个参数的含义如下: b—滤波器系数。对于一个 n阶的 FIR滤波器,其 n+1个滤波器系数可表示为: b(z)=b(1)+b(2)z-1+…+b(n+1)z-n。 n—滤波器阶数。 Wn—截止频率,0≤Wn≤1,Wn=1对应于采样频率的一半。当设计带通和带阻滤波器时,Wn=[W1 W2],W1≤ω≤W2。 ftype—当指定 ftype时,可设计高通和带阻滤波器。Ftype=high时,设计高通 FIR滤波 器;ftype=stop时设计带阻 FIR滤波器。低通和带通 FIR滤波器无需输入 ftype参数。 Window—窗函数。窗函数的长度应等于 FIR滤波器系数个数,即阶数n+1。