4.带阻滤波器:
Wp和Ws为二元矢量且Wp>Ws,如Wp=[0.1,0.8],Ws=[0.2,0.7]。
函数三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=[W1W2],W1≤ω≤W2。
ftype—当指定ftype时,可设计高通和带阻滤波器。
Ftype=high时,设计高通FIR滤波器;ftype=stop时设计带阻FIR滤波器。
低通和带通FIR滤波器无需输入ftype参数。
Window—窗函数。
窗函数的长度应等于FIR滤波器系数个数,即阶数n+1。
函数4窗函数的使用:
在MATLAB下,这些窗函数分别为:
1.矩形窗:
w=boxcar(n),产生一个n点的矩形窗函数。
2.三角窗:
w=triang(n),产生一个n点的三角窗函数。
当n为奇数时,三角窗系数为w(k)=
当n为偶数时,三角窗系数为w(k)=
3.巴特利特窗:
w=Bartlett(n),产生一个n点的巴特利特窗函数。
巴特利特窗系数为w(k)=
巴特利特窗与三角窗非常相似。
巴特利特窗在取样点1和n上总以零结束,而三角窗在这些点上并不为零。
实际上,当n为奇数时bartlett(n)的中心n-2个点等效于triang(n-2)。
4.汉明窗:
w=hamming(n),产生一个n点的汉明窗函数。
汉明窗系数为w(k+1)=0.54-0.46cos()k=0,…,n-1
5.汉宁窗:
w=hanning(n),产生一个n点的汉宁窗函数。
汉宁窗系数为w(k)=0.5[1-cos()]k=1,…,n
6.布莱克曼窗:
w=Blackman(n),产生一个n点的布莱克曼窗函数。
布莱克曼窗系数为w(k)=0.42-0.5cos(2π)+0.8cos(4π)]k=1,…,n
与等长度的汉明窗和汉宁窗相比,布莱克曼窗的主瓣稍宽,旁瓣稍低。
7.凯泽窗:
w=Kaiser(n,beta),产生一个n点的凯泽窗数,其中beta为影响窗函数旁瓣的β参数,其最小的旁瓣抑制α与β的关系为:
0.1102(α-0.87)α>50
β=0.5842(α-21)0.4+0.07886(α-21)21≤α≤50
0α<21
增加β可使主瓣变宽,旁瓣的幅度降低。
8.契比雪夫窗:
w=chebwin(n,r)产生一个n点的契比雪夫窗函数。
其傅里叶变换后的旁瓣波纹低于主瓣r个db数。
Matlab中关于FIR滤波器设计的命令收藏
一、产生窗函数的文件有八个:
1.bartlett(三角窗);——两端为零
2.blackman(布莱克曼窗);
3.boxcar(矩形窗);
4.hamming(哈明窗);
5.hanning(汉宁窗);
6.triang(三角窗);——两端不为零
7.chebwin(切比雪夫窗);
8.kaiser(凯赛窗);
二、fir1.m用“窗函数法”设计FIRDF。
调用格式:
(1)b=fir1(N,Wn);
(2)b=fir1(N,Wn,‘high’);
(3)b=fir1(N,Wn,‘stop’);
N:
阶次,滤波器长度为N+1;
Wn:
通带截止频率,其值在0~1之间,1对应Fs/2;
b:
滤波器系数。
格式
(1):
若Wn为标量,则设计低通滤波器,若Wn是1×2的向量,则用来设计带通滤波器,若Wn是1×L的向量,则可用来设计L带滤波器。
这时,格式
(1)要改为:
b=fir1(N,Wn,'DC-1'),或b=fir1(N,Wn,'DC-0')。
前者保证第一个带为通带,后者保证第一个带为阻带。
格式
(2):
用来设计高通滤波器。
格式(3):
用来设计带阻滤波器。
在上述所有格式中,若不指定窗函数的类型,fir1自动选择Hamming窗。
三、fir2.m本文件采用“窗函数法”设计具有任意幅频相应的FIR数字滤波器。
其调用格式是:
b=fir1(N,F,M);
F是频率向量,其值在0~1之间,M是和F相对应的所希望的幅频相应。
如同fir1,缺省时自动选用Hamming窗。
例:
设计一多带滤波器,要求频率在0.2~0.3,0.6~0.8之间为1,其余处为零。
四、remez.m设计Chebyshev最佳一致逼近FIR滤波器、Hilbert变换器和差分器。
其调用格式是:
(1)b=remez(N,F,A);
(2)b=remez(N,F,A,W);
(3)b=remez(N,F,A,W,‘Hilbert’);
(4)b=remez(N,F,A,W,‘'differentiator')
N是给定的滤波器的阶次,b是设计的滤波器的系数,其长度为N+1;
F是频率向量,A是对应F的各频段上的理想幅频响应,W是各频段上的加权向量。
F、A及W的指定方式和例8.4.1和8.4.2所讨论过的一样,唯一的差别是F的范围为0~1,而非0~0.5,1对应抽样频率的一半。
需要指出的是,若b的长度为偶数,设计高通和带阻滤波器时有可能出现错误,因此,最好保证b的长度为奇数,也即N应为偶数。
五、remezord.m本文件用来确定在用Chebyshev最佳一致逼近设计FIR滤波器时所需要的滤波器阶次。
其调用格式是:
[N,Fo,Ao,W]=remezord(F,A,DEV,Fs)。
F、A的含意同文件remez,DEV是通带和阻带上的偏差;
输出的是适合要求的滤波器阶次N、频率向量Fo、幅度向量Ao和加权向量W。
若设计者事先不能确定要设计的滤波器的阶次,那么,调用remezord后,就可利用这一族参数调用remez,即b=remez(N,Fo,Ao,W),从而设计出所需要滤波器。
因此,remez和remezord常结合起来使用。
需要说明的是,remezord给出的阶次N有可能偏低,这时适当增加N即可;另外,最好判断一下,若N为奇数,就令其加一,使其变为偶数,这样b的长度为奇数。
例:
fedge=[8001000];%表示频率向量,用于低通滤波器的通带截止和阻带起始
mval=[10];%对应fedge各频率向量上的理想幅频响应
dev=[0.05590.01];%通带和阻带上的偏差
fs=4000;%抽样频率
[N,fpts,mag,wt]=remezord(fedge,mval,dev,fs);
%由remezord求得滤波器的阶次N、频率向量fpts、幅度向量mag和加权向量wt
b=remez(N,fpts,mag,wt);
[h,w]=freqz(b,1,256);
plot(w*2000/pi,20*log10(abs(h)));grid;
xlabel('频率/Hz')
ylabel('幅度/dB')
六、其他
1、firls.m用最小平方法设计线性相位FIR滤波器,可设计任意给定的理想幅频响应;
2、fircls.m用带约束的最小平方法设计线性相位FIR滤波器,可设计任意给定的理想幅频响应;
3、fircls1.m用带约束的最小平方方法设计线性相位FIR低通和高通滤波器。
4、sgolay.m用来设计Savitzky-GolayFIR平滑滤波器,其原理见9.1.1节
5、firrcos.m用来设计低通线性相位FIR滤波器,其过渡带为余弦函数形状。
一、在MATLAB中的窗函数,十分简单:
(1)矩形窗(RectangleWindow) 调用格式:
w=boxcar(n),根据长度n产生一个矩形窗w。
(2)三角窗(TriangularWindow) 调用格式:
w=triang(n),根据长度n产生一个三角窗w。
(3)汉宁窗(HanningWindow) 调用格式:
w=hanning(n),根据长度n产生一个汉宁窗w。
(4)海明窗(HammingWindow) 调用格式:
w=hamming(n),根据长度n产生一个海明窗w。
(5)布拉克曼窗(BlackmanWindow)调用格式:
w=blackman(n),根据长度n产生一个布拉克曼窗w。
(6)恺撒窗(KaiserWindow) 调用格式:
w=kaiser(n,beta),根据长度n和影响窗函数旁瓣的β参数产生一个恺撒窗w。
二、基于窗函数的FIR滤波器设计
利用MATLAB提供的函数fir1来实现
调用格式:
fir1(n,Wn,’ftype’,Window),n为阶数、Wn是截止频率(如果输入是形如[W1W2]的矢量时,本函数将设计带通滤波器,其通带为W1<ω<>
[例]设计一个长度为8的线性相位FIR滤波器。
其理想幅频特性满足
用矩形窗:
Window=boxcar(8);
b=fir1(7,0.4,Window);
freqz(b,1)
用blackman窗:
Window=blackman(8);
b=fir1(7,0.4,Window);
freqz(b,1)
[例]设计线性相位带通滤波器,其长度N=15,上下边带截止频率分别为W1=0.3π,w2=0.5π
Window=blackman(16);
b=fir1(15,[0.30.5],Window);
freqz(b,1)
linspace是用来生成一个等差数列的lin是linear的缩写
x=linspace(a,b,n)
就是将生成一个以a和b为断点共划分为n个区间的向量,比如
>>linspace(1,10,5)
ans=
1.00003.25005.50007.750010.0000
AWGN:
在某一信号中加入高斯白噪声
y=awgn(x,SNR)在信号x中加入高斯白噪声。
信噪比SNR以dB为单位。
x的强度假定为0dBW。
如果x是复数,就加入复噪声。
y=awgn(x,SNR,SIGPOWER)如果SIGPOWER是数值,则其代表以dBW为单位的信号强度;如果SIGPOWER为'measured',则函数将在加入噪声之前测定信号强度。
y=awgn(x,SNR,SIGPOWER,STATE)重置RANDN的状态。
y=awgn(…,POWERTYPE)指定SNR和SIGPOWER的单位。
POWERTYPE可以是'dB'或'linear'。
最佳FIR滤波器设计-使用remezord,remez
REMEZ和REMEZORD常用来设计最佳滤波器,其中REMZORD用来计算滤波器的阶数。
使用例子如下:
采用频率fs:
8000Hz
通带截至频率:
1500Hz(归一化后为0.375)
阻带截至频率:
2000Hz(归一化后为0.5)
通带波纹:
0.01(即0.1737dB)
阻带衰减:
0.1(即20dB)
[n,fo,mo,w]=remezord([15002000],[10],[0.010.1],8000);
b=remez(n,fo,mo,w);
%[hw]=freqz(b,1,128);
%plot(w/pi,abs(h));
freqz(b,1,128);
得到的输出响应如下图:
FIR滤波器设计(无反馈)
fir2函数:
基于窗函数的FIR滤波器设计——任意频率响应
语法格式:
b=fir2(n,f,m)
b=fir2(n,f,m,window)
b=fir2(n,f,m,npt)
b=fir2(n,f,m,npt,window)
f频率点矢量(属于[01]区间)
n滤波器阶数
m幅度矢量
npt指定加窗的点数
windows指定的窗函数
2,FIR滤波器设计(无反馈)
fir2函数:
基于窗函数的FIR滤波器设计——任意频率响应
例:
设计一个30阶的带通滤波器
n=30;
f=[00.30.71];
m=[0110];
b=fir2(n,f,m);
freqz(b,1,512);
butter函数:
Butterworth滤波器设计
例2:
从多种频率成分叠加的信号中,提取出某一频带的信号.例如:
从500Hz,1kHz,2kHz三种频率成分叠加的信号中,提取出频带为800Hz~1.5kHz的信号.要求:
阻带衰减不低于10dB,通带衰减不高于3dB.
%先选择滤波器阶数
Fs=5000;%Fs≥2*2kHz
Wp=[8001500]/(Fs/2);
Ws=[6001700]/(Fs/2);
Rp=3;Rs=10;
[n,Wn]=buttord(Wp,Ws,Rp,Rs);
%再进行滤波器设计
[b,a]=butter(n,Wn,'bandpass');
z=filter(b,a,y);
%设计低通滤波器:
[N,Wc]=buttord()
%估算得到Butterworth低通滤波器的最小阶数N和3dB截止频率Wc
[a,b]=butter(N,Wc);%设计Butterworth低通滤波器
[h,f]=freqz();%求数字低通滤波器的频率响应
figure
(2);%打开窗口2
subplot(221);%图形显示分割窗口
plot(f,abs(h));%绘制Butterworth低通滤波器的幅频响应图
title('巴氏低通滤波器');
grid;%绘制带网格的图像
sf=filter(a,b,s);%叠加函数S经过低通滤波器以后的新函数
subplot(222);
plot(t,sf);%绘制叠加函数S经过低通滤波器以后的时域图形
xlabel('时间(seconds)');
ylabel('时间按幅度');
SF=fft(sf,256);%对叠加函数S经过低通滤波器以后的新函数进行256点的基—2快速傅立叶变换
%%%%%%%%%%%%%%%%%%%%%%%%%%w=%新信号角频率
subplot(223);
plot();%绘制叠加函数S经过低通滤波器以后的频谱图
title('低通滤波后的频谱图');
%设计高通滤波器
[N,Wc]=buttord()
%估算得到Butterworth高通滤波器的最小阶数N和3dB截止频率Wc
[a,b]=butter(N,Wc,'high');%设计Butterworth高通滤波器
[h,f]=freqz();%求数字高通滤波器的频率响应
figure(3);
subplot(221);
plot(f,abs(h));%绘制Butterworth高通滤波器的幅频响应图
title('巴氏高通滤波器');
grid;%绘制带网格的图像
sf=filter();%叠加函数S经过高通滤波器以后的新函数
subplot(222);
plot(t,sf);;%绘制叠加函数S经过高通滤波器以后的时域图形
xlabel('Time(seconds)');
ylabel('Timewaveform');
w;%新信号角频率
subplot(223);
plot();%绘制叠加函数S经过高通滤波器以后的频谱图
title('高通滤波后的频谱图');
%设计带通滤波器
[N,Wc]=buttord()
%估算得到Butterworth带通滤波器的最小阶数N和3dB截止频率Wc
[a,b]=butter(N,Wc);%设计Butterworth带通滤波器
[h,f]=freqz();%求数字带通滤波器的频率响应
figure(4);
subplot(221);
plot(f,abs(h));%绘制Butterworth带通滤波器的幅频响应图
title('butterbandpassfilter');
grid;%绘制带网格的图像
sf=filter(a,b,s);%叠加函数S经过带通滤波器以后的新函数
subplot(222);
plot(t,sf);%绘制叠加函数S经过带通滤波器以后的时域图形
xlabel('Time(seconds)');
ylabel('Timewaveform');
SF=fft();%对叠加函数S经过带通滤波器以后的新函数进行256点的基—2快速傅立叶变换
%%%%%%%%%%%%%%%%%%w=(%新信号角频率
subplot(223);
plot();%绘制叠加函数S经过带通滤波器以后的频谱图
title('带通滤波后的频谱图');
abs(x):
纯量的绝对值或向量的长度
angle(z):
复数z的相角(Phaseangle)
sqrt(x):
开平方
real(z):
复数z的实部
imag(z):
复数z的虚部
conj(z):
复数z的共轭复数
round(x):
四舍五入至最近整数
fix(x):
无论正负,舍去小数至最近整数
floor(x):
地板函数,即舍去正小数至最近整数
ceil(x):
天花板函数,即加入正小数至最近整数
rat(x):
将实数x化为分数表示
rats(x):
将实数x化为多项分数展开
sign(x):
符号函数(Signumfunction)。
当x<0时,sign(x)=-1;
当x=0时,sign(x)=0;
当x>0时,sign(x)=1。
rem(x,y):
求x除以y的馀数
gcd(x,y):
整数x和y的最大公因数
lcm(x,y):
整数x和y的最小公倍数
exp(x):
自然指数
pow2(x):
2的指数
log(x):
以e为底的对数,即自然对数或
log2(x):
以2为底的对数
log10(x):
以10为底的对数
二、MATLAB常用的三角函数
sin(x):
正弦函数
cos(x):
馀弦函数
tan(x):
正切函数
asin(x):
反正弦函数
acos(x):
反馀弦函数
atan(x):
反正切函数
atan2(x,y):
四象限的反正切函数
sinh(x):
超越正弦函数
cosh(x):
超越馀弦函数
tanh(x):
超越正切函数
asinh(x):
反超越正弦函数
acosh(x):
反超越馀弦函数
atanh(x):
反超越正切函数
三、适用於向量的常用函数有:
min(x):
向量x的元素的最小值
max(x):
向量x的元素的最大值
mean(x):
向量x的元素的平均值
median(x):
向量x的元素的中位数
std(x):
向量x的元素的标准差
diff(x):
向量x的相邻元素的差
sort(x):
对向量x的元素进行排序(Sorting)
length(x):
向量x的元素个数
norm(x):
向量x的欧氏(Euclidean)长度
sum(x):
向量x的元素总和
prod(x):
向量x的元素总乘积
cumsum(x):
向量x的累计元素总和
cumprod(x):
向量x的累计元素总乘积
dot(x,y):
向量x和y的内积
cross(x,y):
向量x和y的外积
四、MATLAB的永久常数
i或j:
基本虚数单位(即)
eps:
系统的浮点(Floating-point)精确度
inf:
无限大,例如1/0
nan或NaN:
非数值(Notanumber),例如0/0
pi:
圆周率p(=3.1415926...)
realmax:
系统所能表示的最大数值
realmin:
系统所能表示的最小数值
nargin:
函数的输入引数个数
nargin:
函数的输出引数个数
五、MATLAB基本绘图函数
plot:
x轴和y轴均为线性刻度(Linearscal