完整word版,各类滤波器的MATLAB程序清单
- 格式:doc
- 大小:575.53 KB
- 文档页数:12
MATLAB小波变换指令及其功能介绍1 一维小波变换的 Matlab 实现(1) dwt函数功能:一维离散小波变换格式:[cA,cD]=dwt(X,'wname')[cA,cD]=dwt(X,Lo_D,Hi_D)别可以实现一维、二维和 N 维 DFT说明:[cA,cD]=dwt(X,’wname’) 使用指定的小波基函数’wname’ 对信号X 进行分解,cA、cD 分别为近似分量和细节分量;[cA,cD]=dwt(X,Lo_D,Hi_D) 使用指定的滤波器组 Lo_D、Hi_D 对信号进行分解。
(2) idwt 函数功能:一维离散小波反变换格式:X=idwt(cA,cD,'wname')X=idwt(cA,cD,Lo_R,Hi_R)X=idwt(cA,cD,'wname',L)函数 fft、fft2 和 fftn 分X=idwt(cA,cD,Lo_R,Hi_R,L)说明:X=idwt(cA,cD,’wname’) 由近似分量 cA 和细节分量 cD 经小波反变换重构原始信号 X 。
'wname’为所选的小波函数X=idwt(cA,cD,Lo_R,Hi_R)用指定的重构滤波器 Lo_R 和 Hi_R 经小波反变换重构原始信号 X 。
X=idwt(cA,cD,’wname',L) 和 X=idwt(cA,cD,Lo_R,Hi_R,L)指定返回信号 X 中心附近的 L 个点。
2 二维小波变换的 Matlab 实现二维小波变换的函数别可以实现一维、二维和 N 维 DFT函数名函数功能—-————---—--—---——---—---—-—---—-——----——-----—————dwt2 二维离散小波变换wavedec2 二维信号的多层小波分解idwt2 二维离散小波反变换waverec2 二维信号的多层小波重构wrcoef2 由多层小波分解重构某一层的分解信号upcoef2 由多层小波分解重构近似分量或细节分量detcoef2 提取二维信号小波分解的细节分量appcoef2 提取二维信号小波分解的近似分量upwlev2 二维小波分解的单层重构dwtpet2 二维周期小波变换idwtper2 二维周期小波反变换—-----—-—-—-—-—-—--—-—-------—-——-—-————-———-—-——-——-—-----(1) wcodemat 函数功能:对数据矩阵进行伪彩色编码函数 fft、fft2 和 fftn 分格式:Y=wcodemat(X,NB,OPT,ABSOL)Y=wcodemat(X,NB,OPT)Y=wcodemat(X,NB)Y=wcodemat(X)说明:Y=wcodemat(X,NB,OPT,ABSOL)返回数据矩阵 X 的编码矩阵 Y ;NB 伪编码的最大值,即编码范围为 0~NB,缺省值 NB=16;OPT 指定了编码的方式(缺省值为 'mat'),即:别可以实现一维、二维和 N 维 DFTOPT='row’ ,按行编码OPT='col' ,按列编码OPT=’mat’ ,按整个矩阵编码函数 fft、fft2 和 fftn 分ABSOL 是函数的控制参数(缺省值为’1'),即:ABSOL=0 时,返回编码矩阵ABSOL=1 时,返回数据矩阵的绝对值 ABS(X)1. 离散傅立叶变换的Matlab实现(2) dwt2 函数功能:二维离散小波变换格式:[cA,cH,cV,cD]=dwt2(X,'wname')[cA,cH,cV,cD]=dwt2(X,Lo_D,Hi_D)说明:[cA,cH,cV,cD]=dwt2(X,'wname’)使用指定的小波基函数'wname’ 对二维信号 X 进行二维离散小波变幻;cA,cH,cV,cD 分别为近似分量、水平细节分量、垂直细节分量和对角细节分量;[cA,cH,cV,cD]=dwt2(X,Lo_D,Hi_D)使用指定的分解低通和高通滤波器 Lo_D 和 Hi_D 分解信号 X 。
matlab自带的滤波器函数
Matlab自带的滤波器函数可以用于对信号进行滤波处理,常用的函数有:
1. fir1函数:设计一阶低通、高通、带通、带阻滤波器的FIR 数字滤波器,可自定义通带和阻带的截止频率。
2. cheby1函数:设计ChebyshevI型低通、高通、带通、带阻数字滤波器,可自定义通带和阻带的截止频率和最大通带波纹。
3. butter函数:设计Butterworth型低通、高通、带通、带阻数字滤波器,可自定义通带和阻带的截止频率和滤波器阶数。
4. filtfilt函数:对信号进行双向滤波处理,可避免滤波后信号的相位畸变和滞后。
这些函数可以在Matlab的Signal Processing Toolbox中找到,可根据需要选择合适的函数进行滤波处理。
- 1 -。
各类滤波器的MATLAB程序一、理想低通滤波器IA=imread('lena.bmp');[f1,f2]=freqspace(size(IA),'meshgrid');Hd=ones(size(IA));r=sqrt(f1.^2+f2.^2);Hd(r>0.2)=0;Y=fft2(double(IA));Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA));subplot(2,2,2),imshow(uint8(Ia));figuresurf(Hd,'Facecolor','interp','Edgecolor','none','Facelighting','phong');二、理想高通滤波器IA=imread('lena.bmp');[f1,f2]=freqspace(size(IA),'meshgrid');Hd=ones(size(IA));r=sqrt(f1.^2+f2.^2);Hd(r<0.2)=0;Y=fft2(double(IA));Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya));figuresubplot(2,2,1),imshow(uint8(IA));subplot(2,2,2),imshow(uint8(Ia));figuresurf(Hd,'Facecolor','interp','Edgecolor','none','Facelighting','phong');三、B utterworth低通滤波器IA=imread('lena.bmp');[f1,f2]=freqspace(size(IA),'meshgrid');D=0.3;r=f1.^2+f2.^2;n=4;for i=1:size(IA,1)for j=1:size(IA,2)t=r(i,j)/(D*D);Hd(i,j)=1/(t^n+1);endendY=fft2(double(IA));Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya));figuresubplot(2,2,1),imshow(uint8(IA));subplot(2,2,2),imshow(uint8(Ia));figuresurf(Hd,'Facecolor','interp','Edgecolor','none','Facelighting','phong');四、B utterworth高通滤波器IA=imread('lena.bmp');[f1,f2]=freqspace(size(IA),'meshgrid');D=0.3;r=f1.^2+f2.^2;n=4;for i=1:size(IA,1)for j=1:size(IA,2)t=(D*D)/r(i,j);Hd(i,j)=1/(t^n+1);endendY=fft2(double(IA));Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya));figuresubplot(2,2,1),imshow(uint8(IA));subplot(2,2,2),imshow(uint8(Ia));figuresurf(Hd,'Facecolor','interp','Edgecolor','none','Facelighting','phong');五、高斯低通滤波器IA=imread('lena.bmp');IB=imread('babarra.bmp');[f1,f2]=freqspace(size(IA),'meshgrid');D=100/size(IA,1);r=f1.^2+f2.^2;Hd=ones(size(IA));for i=1:size(IA,1)for j=1:size(IA,2)t=r(i,j)/(D*D);Hd(i,j)=exp(-t);endendY=fft2(double(IA));Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya));figuresubplot(2,2,1),imshow(uint8(IA));subplot(2,2,2),imshow(uint8(Ia));figuresurf(Hd,'Facecolor','interp','Edgecolor','none','Facelighting','phong');六、高斯高通滤波器IA=imread('lena.bmp');IB=imread('babarra.bmp');[f1,f2]=freqspace(size(IA),'meshgrid');%D=100/size(IA,1);D=0.3;r=f1.^2+f2.^2;for i=1:size(IA,1)for j=1:size(IA,2)t=r(i,j)/(D*D);Hd(i,j)=1-exp(-t);endendY=fft2(double(IA));Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya));figuresubplot(2,2,1),imshow(uint8(IA));subplot(2,2,2),imshow(uint8(Ia));figuresurf(Hd,'Facecolor','interp','Edgecolor','none','Facelighting','phong');七、梯形低通滤波器IA=imread('lena.bmp');IB=imread('babarra.bmp');[f1,f2]=freqspace(size(IA),'meshgrid');%D=100/size(IA,1);D0=0.1;D1=0.4;r=sqrt(f1.^2+f2.^2);Hd=zeros(size(IA));Hd(r<D0)=1;for i=1:size(IA,1)for j=1:size(IA,2)if r(i,j)>=D0 & r(i,j)<=D1Hd(i,j)=(D1-r(i,j))/(D1-D0);endendendY=fft2(double(IA));Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya));figuresubplot(2,2,1),imshow(uint8(IA));subplot(2,2,2),imshow(uint8(Ia));figuresurf(Hd,'Facecolor','interp','Edgecolor','none','Facelighting','phong');八、梯形高通滤波器IA=imread('lena.bmp');IB=imread('babarra.bmp');[f1,f2]=freqspace(size(IA),'meshgrid');%D=100/size(IA,1);D0=0.1;D1=0.4;r=sqrt(f1.^2+f2.^2);Hd=ones(size(IA));Hd(r<D1)=0;for i=1:size(IA,1)for j=1:size(IA,2)if r(i,j)>=D0 & r(i,j)<=D1Hd(i,j)=(D0-r(i,j))/(D0-D1);endendendY=fft2(double(IA));Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya));figuresubplot(2,2,1),imshow(uint8(IA));subplot(2,2,2),imshow(uint8(Ia));figuresurf(Hd,'Facecolor','interp','Edgecolor','none','Facelighting','phong');九、用其他方法编写的理想低通、理想高通、Butterworth低通、同态滤波程序1、理想低通i1=imread('lena.bmp');i2=imnoise(i1,'salt & pepper',0.1);f=double(i2);k=fft2(f);g=fftshift(k);[N1,N2]=size(g);d0=50;u0=floor(N1/2)+1;v0=floor(N2/2)+1;for i=1:N1for j=1:N2d=sqrt((i-u0)^2+(j-v0)^2);if d<=d0h=1;elseh=0;endy(i,j)=g(i,j)*h;endendy=ifftshift(y);E1=ifft2(y);E2=real(E1);figuresubplot(2,2,1),imshow(uint8(i1));subplot(2,2,2),imshow(uint8(i2));subplot(2,2,3),imshow(uint8(E2));2、理想高通i1=imread('lena.bmp');i2=imnoise(i1,'salt & pepper',0.1); f=double(i2);k=fft2(f);g=fftshift(k);[N1,N2]=size(g);n=2;d0=10;u0=floor(N1/2)+1;v0=floor(N2/2)+1;for i=1:N1for j=1:N2d=sqrt((i-u0)^2+(j-v0)^2);if d<=d0h=0;else h=1;endy(i,j)=g(i,j)*h;endendy=ifftshift(y);E1=ifft2(y);E2=real(E1);figuresubplot(2,2,1),imshow(uint8(i1)); subplot(2,2,2),imshow(uint8(i2)); subplot(2,2,3),imshow(uint8(E2));3、Butterworth低通i1=imread('lena.bmp');i2=imnoise(i1,'salt & pepper',0.1); f=double(i2);k=fft2(f);g=fftshift(k);[N1,N2]=size(g);n=2;d0=50;u0=floor(N1/2)+1;v0=floor(N2/2)+1;for i=1:N1for j=1:N2d=sqrt((i-u0)^2+(j-v0)^2);h=1/(1+(d/d0)^(2*n));y(i,j)=g(i,j)*h;endendy=ifftshift(y);E1=ifft2(y);E2=real(E1);figuresubplot(2,2,1),imshow(uint8(i1)); subplot(2,2,2),imshow(uint8(i2)); subplot(2,2,3),imshow(uint8(E2));4、同态滤波I=rgb2gray(imread('fabric00.bmp')); [M,N]=size(I);T=double(I);L=log(T);F=fft2(L);A=2;B=0.3;for i=1:Mfor j=1:ND(i,j)=((i-M/2)^2+(j-N/2)^2);endendc=1.1;%锐化参数D0=max(M,N);H=(A-B)*(1-exp(c*(-D/(D0^2))))+B;F=F.*H;F=ifft2(F);Y=exp(F);figuresubplot(1,2,1),imshow(I);subplot(1,2,2),imshow(uint8(real(Y)));十、十一、Gabor滤波器。
matlab滤波器程序wc=(400/1000)*pi;%求截止频率w1=boxcar(81);%窗函数w2=triang(81);w3=hamming(81);w4=hanning(81);w5=bartlett(81);w6=blackman(81);w7=chebwin(81,30);w8=kaiser(81,7.856);n=1:1:81;hd=sin(wc*(n-41))./(pi*(n-41)); %求h(d)hd(41)=wc/pi;h1=hd.*w1';%加窗h2=hd.*w2';h3=hd.*w3';h4=hd.*w4';h5=hd.*w5';h6=hd.*w6';h7=hd.*w7';h8=hd.*w8';[mag1,rad]=freqz(h1);%求幅频特性曲线[mag2,rad]=freqz(h2);[mag3,rad]=freqz(h3);[mag4,rad]=freqz(h4);[mag5,rad]=freqz(h5);[mag6,rad]=freqz(h6);[mag7,rad]=freqz(h7);[mag8,rad]=freqz(h8);figure(1);%画幅频特性曲线plot(rad,20*log10(abs(mag1)));xlabel('Normalized Frequency(rad)');ylabel('Normaliaed Magnitude(dB)');axis([0,3,-80,0]);title('利用矩形窗设计的数字滤波器');grid on;figure(2);plot(rad,20*log10(abs(mag2)));xlabel('Normalized Frequency(rad)');ylabel('Normaliaed Magnitude(dB)');axis([0,3,-80,0]);title('利用三角窗设计的数字滤波器');grid on;figure(3);plot(rad,20*log10(abs(mag3)));xlabel('Normalized Frequency(rad)');ylabel('Normaliaed Magnitude(dB)');axis([0,3,-80,0]);title('利用海明设计的数字滤波器');grid on;figure(4);plot(rad,20*log10(abs(mag4)));xlabel('Normalized Frequency(rad)');ylabel('Normaliaed Magnitude(dB)');axis([0,3,-80,0]);title('利用汉宁窗设计的数字滤波器');grid on;figure(5);plot(rad,20*log10(abs(mag5)));xlabel('Normalized Frequency(rad)');ylabel('Normaliaed Magnitude(dB)');axis([0,3,-80,0]);title('利用巴特里特窗设计的数字滤波器');grid on;figure(6);plot(rad,20*log10(abs(mag6)));xlabel('Normalized Frequency(rad)');ylabel('Normaliaed Magnitude(dB)');axis([0,3,-100,0]);title('利用布拉克曼窗设计的数字滤波器');grid on;figure(7);plot(rad,20*log10(abs(mag7)));xlabel('Normalized Frequency(rad)');ylabel('Normaliaed Magnitude(dB)');axis([0,3,-100,0]);title('利用切比雪夫窗设计的数字滤波器');grid on;figure(8);plot(rad,20*log10(abs(mag8)));xlabel('Normalized Frequency(rad)');ylabel('Normaliaed Magnitude(dB)');axis([0,3,-100,0]);title('利用凯塞窗设计的数字滤波器');grid on;MATLAB参考程序和仿真内容%*******************************************************************% %mode: 1--巴特沃斯低通;2--巴特沃斯高通;3--巴特沃斯带通;4--巴特沃斯带阻% 5--契比雪夫低通;6--契比雪夫高通;7--契比雪夫带通;8--契比雪夫带阻%fp1,fp2:通带截止频率,当高通或低通时只有fp1有效%fs1, fs2:阻带截止频率,当高通或低通时只有fs1有效%rp: 通带波纹系数%as: 阻带衰减系数%sample: 采样率%h: 返回设计好的滤波器系数%*******************************************************************% function[b,a]=iirfilt(mode,fp1,fp2,fs1,fs2,rp,as,sample)wp1=2*fp1/sample;wp2=2*fp2/sample;ws1=2*fs1/sample;ws2=2*fs2/sample;%得到巴特沃斯滤波器的最小阶数N和3bd频率wnif mode<3[N,wn]=buttord(wp1,ws1,rp,as);elseif mode<5[N,wn]=buttord([wp1 wp2],[ws1 ws2],rp,as);%得到契比雪夫滤波器的最小阶数N和3bd频率wnelseif mode<7[N,wn]=cheb1ord(wp1,ws1,rp,as);else[N,wn]=cheblord([wp1 wp2],[ws1 ws2],rp,as);end%得到滤波器系数的分子b和分母aif mode= =1[b,a]=butter(N,wn);endif mode= =2[b,a]=butter(N,wn,/high/);endif mode= =3[b,a]=butter(N,wn);endif mode= =4[b,a]=butter(N,wn,/stop/);endif mode= =5[b,a]=cheby1(N,rp,wn);endif mode= =6[b,a]=cheby1(N,rp,wn,/high/);endif mode= =7[b,a]=cheby1(N,rp,wn);endif mode= =8[b,a]=cheby1(N,rp,wn,/stop/);endset(gcf,/menubar/,menubar);freq_response=freqz(b,a);magnitude=20*log10(abs(freq_response));m=0:511;f=m*sample/(2*511);subplot(3,1,1);plot(f,magnitude);grid; %幅频特性axis([0 sample/2 1.1*min(magnitude) 1.1*max(magnitude)]);ylabel('Magnitude');xlabel('Frequency-->');phase=angle(freq_response);subplot(3,1,2);plot(f,phase);grid; %相频特性axis([0 sample/2 1.1*min(phase) 1.1*max(phase)]);ylabel('Phase');xlabel('Frequency-->');h=impz(b,a,32); %32点的单位函数响应t=1:32;subplot(3,1,3);stem(t,h);grid;axis([0 32 1.2*min(h) 1.1*max(h)]);ylabel('h(n)');xlabel('n-->');基于MATLAB信号处理工具箱的数字滤波器设计与仿真摘要:传统的数字滤波器的设计过程复杂,计算工作量大,滤波特性调整困难,影响了它的应用。
第二章 MATLAB 语言及应用实验项目实验一 MATLAB 数值计算三、实验内容与步骤1.创建矩阵⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=987654321a(1(2)用(3)用(42.矩阵的运算(1)利用矩阵除法解线性方程组。
⎪⎪⎩⎪⎪⎨⎧=+++=-+-=+++=+-12224732258232432143214321421x x x x x x x x x x x x x x x 将方程表示为AX=B ,计算X=A\B 。
(2)利用矩阵的基本运算求解矩阵方程。
已知矩阵A 和B 满足关系式A -1BA=6A+BA ,计算矩阵B 。
其中⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=7/10004/10003/1A ,Ps: format rata=[1/3 0 0;0 1/4 0;0 0 1/7];b=inv(a)*inv(inv(a)-eye(3))*6*a(3)计算矩阵的特征值和特征向量。
已知矩阵⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--=1104152021X ,计算其特征值和特征向量。
(4)Page:322利用数学函数进行矩阵运算。
已知传递函数G(s)=1/(2s+1),计算幅频特性Lw=-20lg(1)2(2w )和相频特性Fw=-arctan(2w),w 的范围为[0.01,10],按对数均匀分布。
3.多项式的运算(1)多项式的运算。
已知表达式G(x)=(x-4)(x+5)(x 2-6x+9),展开多项式形式,并计算当x 在[0,20]内变化时G(x)的值,计算出G(x)=0的根。
Page 324(2)多项式的拟合与插值。
将多项式G(x)=x 4-5x 3-17x 2+129x-180,当x 在[0,20]多项式的值上下加上随机数的偏差构成y1,对y1进行拟合。
对G(x)和y1分别进行插值,计算在5.5处的值。
Page 325 四、思考练习题1.使用logspace 函数创建0~4π的行向量,有20个元素,查看其元素分布情况。
Ps: logspace(log10(0),log10(4*pi),20) (2) sort(c,2) %顺序排列 3.1多项式1)f(x)=2x 2+3x+5x+8用向量表示该多项式,并计算f(10)值. 2)根据多项式的根[-0.5 -3+4i -3-4i]创建多项式。
fs=15000;T= 1/fs;rp=1;rs=40;wp1=0.11*pi;wp2=0.81*pi;ws1=0.31*pi;ws2=0.61*pi;%数字带阻滤波器技术指标wc1=(2/T)*tan(wp1/2);%频率预畸变wc2=(2/T)*tan(wp2/2);wr1=(2/T)*tan(ws1/2);wr2=(2/T)*tan(ws2/2);w0=sqrt(wc1*wc2);B=wc2-wc1;wp=1;%归一化通带截止频率ws=wp*(wr1*B) / (w0^2-wr1^2) ; %归一化阻带截止频率[N,wc]=buttord(wp,ws,rp,rs,'s')%求滤波器阶数和3dB截止频率[Z,P,K]=buttap(N)%设计模拟低通滤波器[Md,Nd]=zp2tf(Z,P,K)%将零极点形式转换为传输函数形式[M,N]=lp2bs(Md,Nd,w0,B)%对低通滤波器进行频率变换,转换为带阻滤波器[h,w]=freqs(M,N);%模拟带阻滤波器的幅频响应plot(w/(2*pi),abs(h));grid;xlabel('频率Hz');ylabel('幅度');title('模拟带阻滤波器');[b,a]=bilinear(M,N,15000)%对模拟滤波器双线性变换figure(1);freqz(b,a);[H,W]=freqz(b,a); %绘出频率响应;axis([0,1,-100,20]);figure(2);plot(W*fs/(2*pi),abs(H));grid on;xlabel('频率/Hz');ylabel('幅值');n=0:199;t=n/fs;x=sin(2*pi*400*t)+3*sin(2*pi*3000*t)+2*sin(2*pi*5000*t);figure(3);subplot(311);plot(t,x);axis([0,0.01,-5,5]);title('输入信号');grid on;y=filter(b,a,x);subplot(312);stem(y,'.');title('输出序列');grid on;ya=y*sinc(fs*(ones(length(n),1)*t-(n/fs)'*ones(1,length(t))));subplot(313);plot(t,ya);axis([0,0.01,-3,3]);title('输出波形');grid on;t=(0:100)/fs;figure(4)fs=1.5*10000;n=(0:100)/fs;f=sin(2*pi*400*t)+3*sin(2*pi*3000*t)+2*sin(2*pi*5000*t);y=fftfilt(b,x);[H1,f1]=freqz(f,[1]);[H2,f2]=freqz(y,[1]);f1=f1/pi*fs/2;f2=f2/pi*fs/2;subplot(2,1,1);plot(f1,abs(H1));title('输入信号的频谱');subplot(2,1,2);plot(f2,abs(H2));title('输出信号的频谱');基于Matlab 的带阻滤波器设计.10.20.30.40.50.60.70.80.91-800-600-400-2000N o r m a l i z e d Fre q u⨯π r a d /s a m p l e Ph a se(d e g r e e s )00.10.20.30.40.50.60.70.80.91-100-50N o r m a l i z e d Fr e q u⨯π r a d /s a m p l e M a g n i tu d e1000200030004000500060007000800000.20.40.60.811.21.4频率/Hz幅值00.0010.0020.0030.0040.0050.0060.0070.0080.0090.01-505输入信号020406080100120140160180200-22输出序列0.0010.0020.0030.0040.0050.0060.0070.0080.0090.01-202输出波形01000200030004000500060007000800050100150200输入信号的频谱010002000300040005000600070008000102030输出信号的频谱N =4wc =1.7947b =0.0186 -0.0410 0.1082 -0.1355 0.1810 -0.1355 0.1082 -0.0410 0.0186a =1.0000 -0.6707 -1.3750 0.5678 1.1964 -0.2996 -0.4631 0.0496 0.0762>。
>> n=31;%定义滤波器阶数32fs=12.8*10^3;fc1=49;fc2=51;w1=2*pi*fc1/fs;w2=2*pi*fc2/fs;%参数转换,将模拟滤波器的技术指标转换为数字滤波器的技术指标window=hanning(n+1);%使用hanning窗函数q=fir1(n,[w1/pi w2/pi],hanning(n+1));%滤波器时域函数,使用标准响应的加窗设计函数fir1w=linspace(0,pi,512);h1=freqz(q,1,512);%进行512个点的傅里叶变换figure(2);plot(w/pi,20*log10(abs(h1)));title('滤波器频谱图');xlabel('频率');ylabel('幅度');grid ;设计FIR低通滤波器,系统频率为50MHz,通带截止频率Fpass为1MHz,阻带截止频率Fstop为4MHz,通带最大衰减Apass为1dB,阻带最小衰减Astop 为30dB。
程序和必要的程序注释谢谢最佳答案只要用一个公式就行。
library IEEE;use IEEE.STD_LOGIC_1164.ALL;use IEEE.STD_LOGIC_ARITH.ALL;use IEEE.STD_LOGIC_UNSIGNED.ALL;entity fir isPort (clk: in std_logic;reset: in std_logic;inpx: in std_logic_vector(11 downto 0);outy: out std_logic_vector(11 downto 0));end fir;architecture beh of fir issignal x0,x1,x2,x3: std_logic_vector(11 downto 0);constant c0:integer :=-1234*32768/1000;constant c1:integer :=2345*32768/10000;constant c2:integer :=5*32768;constant c3:integer :=-3*32768/10000;signal p0,p1,p2,p3:integer;signal sum: integer;beginsample_delay_line:process(clk)beginif rising_edge(clk) thenif reset='1' thenx3 <=(others=>'0');x2 <=(others=>'0');x1 <=(others=>'0');x0 <=(others=>'0');elsex3 <=x2;x2 <=x1;x1 <=x0;x0 <=inpx;end if;end if;end process;p0 <= conv_integer(x0)*c0;p1 <= conv_integer(x1)*c1;p2 <= conv_integer(x2)*c2;p3 <= conv_integer(x3)*c3;sum <=p0+p1+p2+p3;outy <=conv_std_logic_vector(sum/32768,12);end beh;4.1 数字滤波器简介数字滤波器是一种用来过滤时间离散信号的数字系统,通过对抽样数据进行数学处理来达到频域滤波的目的。
滤波器设计MATLAB滤波器的设计在信号处理中具有重要的作用,可以用于去除噪声、增强信号等。
MATLAB是一种强大的工具,可以用于滤波器设计和分析。
本文将介绍如何使用MATLAB进行滤波器设计,并通过示例展示具体的过程。
在MATLAB中,可以使用信号处理工具箱提供的函数来设计滤波器。
常用的函数有:- `fir1`:设计FIR滤波器。
- `butter`:设计巴特沃斯滤波器。
- `cheby1`:设计切比雪夫I型滤波器。
- `cheby2`:设计切比雪夫II型滤波器。
- `ellip`:设计椭圆滤波器。
这些函数的输入参数包括滤波器类型、阶数、截止频率等。
根据具体的需求选择不同的函数来设计滤波器。
下面以设计一个低通滤波器为例,演示如何使用MATLAB进行滤波器设计。
首先,创建一个信号作为输入。
可以使用`sin`函数生成一个正弦信号作为示例。
代码如下:```matlabfs = 1000; % 采样率t = 0:1/fs:1; % 时间向量f=50;%信号频率x = sin(2*pi*f*t); % 输入信号```接下来,使用`fir1`函数设计一个低通滤波器。
该函数的输入参数`n`表示滤波器的阶数,`Wn`表示归一化的截止频率。
代码如下:```matlabn=50;%滤波器阶数Wn=0.2;%截止频率b = fir1(n, Wn);```然后,使用`filter`函数对输入信号进行滤波。
该函数的输入参数是滤波器的系数和输入信号。
代码如下:```matlaby = filter(b, 1, x);```最后,绘制原始信号和滤波后的信号的时域和频域波形。
代码如下:```matlab%时域波形subplot(2, 1, 1)plot(t, x)hold onplot(t, y)legend('原始信号', '滤波后信号') xlabel('时间 (s)')ylabel('幅值')title('时域波形')%频域波形subplot(2, 1, 2)f = linspace(-fs/2, fs/2, length(x)); X = abs(fftshift(fft(x)));Y = abs(fftshift(fft(y)));plot(f, X)hold onplot(f, Y, 'r')legend('原始信号', '滤波后信号') xlabel('频率 (Hz)')ylabel('幅值')title('频域波形')```运行以上代码,可以得到原始信号和滤波后信号的时域和频域波形图。
汉明窗的FIR低通滤波:Fs=22050;[x,FS,bits]=wavread('G:\Users\DP\Desktop\SoundTest.wav');%G:\Users\DP\Desktop \SoundTest.wav G:\Users\DP\Desktop\好铃网-湖人掘金宣传片.wavx=x(:,1);figure(1);subplot(2,1,1);plot(x);%sound(x,FS,bits); %回放语音title('语音信号时域波形图')y=fft(x,3260);f=(FS/3260)*[1:1630];subplot(2,1,2);plot(f(1:1630),abs(y(1:1630)));title('语音信号频谱图');%产生噪声信号并加到语音信号t=0:length(x)-1;zs0=0.05*cos(2*pi*100*t/22050);figure(2);subplot(2,1,1)plot(zs0)title('噪声信号波形');zs1=fft(zs0,1200);%sound(zs0,FS,bits);%回放噪音f=(FS/1200)*[1:600];subplot(2,1,2)plot(f(1:600),abs(zs1(1:600)));title('噪声信号频谱');%sound(x1,FS,bits); %回放加入噪声后的语音y1=fft(x1,1200);figure(3);subplot(2,1,1);plot(x1);f=(FS/1200)*[1:600];subplot(2,1,2);plot(f(1:600),abs(y1(1:600)));title('加入噪声后的信号波形');title('加入噪声后的信号频谱');%加窗滤波wp=0.25*pi;ws=0.3*pi;wdelta=ws-wp;N=ceil(6.6*pi/wdelta); %取整t=0:(size(x1)-1);wn=(0.2+0.3)*pi/2;b=fir1(N,wn/pi,hamming(N+1)); %选择窗函数,并归一化截止频率f1=fftfilt(b,x1);figure(4)freqz(b,1,512)[h1,w1]=freqz(b,1);plot(w1*FS/(2*pi),20*log10(abs(h1)));figure(5)subplot(2,1,1)plot(t,x1)title('滤波前的时域波形');subplot(2,1,2)title('滤波后的时域波形');sound(f1); %播放滤波后的语音信号F0=fft(f1,1024);f=FS*(0:511)/1024;figure(6)y2=fft(x1,1024);subplot(2,1,1);plot(f,abs(y2(1:512))); %画出滤波前的频谱图title('滤波前的频谱')xlabel('Hz');ylabel('fuzhi');subplot(2,1,2)F1=plot(f,abs(F0(1:512))); %画出滤波后的频谱图title('滤波后的频谱')xlabel('Hz');ylabel('fuzhi');巴特沃兹低通滤波器:Fs=22050;[x,FS,bits]=wavread('G:\Users\DP\Desktop\SoundTest.wav');%G:\Users\DP\Desktop \SoundTest.wav G:\Users\DP\Desktop\好铃网-湖人掘金宣传片.wavx=x(:,1);figure(1);subplot(2,1,1);plot(x);%sound(x,FS,bits); %回放语音title('语音信号时域波形图')y=fft(x,3260);f=(FS/1630)*[1:1630];subplot(2,1,2);plot(f(1:1630),abs(y(1:1630)));title('语音信号频谱图');%产生噪声信号并加到语音信号t=0:length(x)-1;zs=0.05*cos(2*pi*10000*t/22050);zs0=0.05*cos(2*pi*10000*t/22050000); figure(2);subplot(2,1,1)plot(zs0)title('噪声信号波形');zs1=fft(zs,1200);%sound(zs,FS,bits); %回放噪音subplot(2,1,2)plot(f(1:600),abs(zs1(1:600)));title('噪声信号频谱');x1=x+zs';sound(x1,FS,bits); %回放加入噪声后的语音y1=fft(x1,1200);figure(3);subplot(2,1,1);plot(x1);title('加入噪声后的信号波形');subplot(2,1,2);plot(f(1:600),abs(y1(1:600)));title('加入噪声后的信号频谱');%低通滤波fp=3000;fs=3500;Fs=22050;rp=1;rs=10;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;Fs1=1;wap=2*tan(wp/2);was=2*tan(ws/2);[N,wc]=buttord(wap,was,rp,rs,'s');[B,A]=butter(N,wc,'s');[Bz,Az]=bilinear(B,A,Fs1);figure(4);[h,w]=freqz(Bz,Az,512,Fs1*22050);plot(w,abs(h));title('巴特沃斯低通滤波器');xlabel('频率(HZ)');ylabel('耗损(dB)');grid on;yd=filter(Bz,Az,x1);figure(5);subplot(2,1,1);plot(yd); title('滤波后信号波形');ydd=fft(yd,1200);subplot(2,1,2);plot(f(1:600),abs(ydd(1:600)));title('滤波后信号频谱');%sound(yd,FS,bits);巴特沃兹高通滤波器:Fs=22050;[x,FS,bits]=wavread('G:\Users\DP\Desktop\SoundTest.wav'); x=x(:,1);figure(1);subplot(2,1,1);plot(x);%sound(x,FS,bits); %回放语音title('语音信号时域波形图')y=fft(x,3260);f=(FS/1630)*[1:1630];subplot(2,1,2);plot(f(1:1630),abs(y(1:1630)));title('语音信号频谱图');%产生噪声信号并加到语音信号t=0:length(x)-1;zs0=0.25*cos(2*pi*100*t/22050);figure(2);subplot(2,1,1)plot(zs0)title('噪声信号波形');zs1=fft(zs0,1200);%sound(zs0,FS,bits); %回放噪音subplot(2,1,2)plot(f(1:600),abs(zs1(1:600)));title('噪声信号频谱');x1=x+zs0';%sound(x1,FS,bits); %回放加入噪声后的语音y1=fft(x1,1200);figure(3);subplot(2,1,1);plot(x1);title('加入噪声后的信号波形');subplot(2,1,2);plot(f(1:600),abs(y1(1:600)));title('加入噪声后的信号频谱');%高通滤波fp=400;fs=300;Fs=22050;rp=1;rs=10;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;T=1;Fs1=1;wap=2*tan(wp/2);was=2*tan(ws/2);[N,wc]=buttord(wap,was,rp,rs,'s');[B,A]=butter(N,wc,'high','s');[Bz,Az]=bilinear(B,A,Fs1);figure(4);[h,w]=freqz(Bz,Az,512,Fs1*22050);plot(w,abs(h));title('巴特沃斯高通滤波器');xlabel('频率(HZ)');ylabel('耗损(dB)');grid on;yd=filter(Bz,Az,x1);figure(5);subplot(2,1,1);plot(yd); title('滤波后信号波形');ydd=fft(yd,1200);subplot(2,1,2);plot(f(1:600),abs(ydd(1:600)));title('滤波后信号频谱');%sound(yd,FS,bits)巴特沃兹带通滤波器:Fs=22050;[x,FS,bits]=wavread('G:\Users\DP\Desktop\SoundTest.wav');x=x(:,1);figure(1);subplot(2,1,1);plot(x);%sound(x,FS,bits); %回放语音title('语音信号时域波形图')y=fft(x,3260);f=(FS/1630)*[1:1630];subplot(2,1,2);plot(f(1:1630),abs(y(1:1630)));title('语音信号频谱图');%产生噪声信号并加到语音信号t=0:length(x)-1;zs0=0.05*cos(2*pi*100*t/22050);figure(2);subplot(2,1,1)plot(zs0)title('噪声信号波形');zs1=fft(zs0,1200);%sound(zs,FS,bits); %回放噪音subplot(2,1,2)plot(f(1:600),abs(zs1(1:600)));title('噪声信号频谱');x1=x+zs0';%sound(x1,FS,bits); %回放加入噪声后的语音y1=fft(x1,1200);figure(3);subplot(2,1,1);plot(x1);title('加入噪声后的信号波形');subplot(2,1,2);plot(f(1:600),abs(y1(1:600)));title('加入噪声后的信号频谱');%带通滤波fp=[800,9000];fs=[700,10000];Fs=22050;rp=1;rs=10;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;T=1;Fs1=1;wap=2*tan(wp/2);was=2*tan(ws/2);[N,wc]=buttord(wap,was,rp,rs,'s');[B,A]=butter(N,wc,'s');[Bz,Az]=bilinear(B,A,Fs1);figure(4);[h,w]=freqz(Bz,Az,512,Fs1*22050);plot(w,abs(h));title('巴特沃斯带通滤波器');xlabel('频率(HZ)');ylabel('耗损(dB)'); grid on;yd=filter(Bz,Az,x1);figure(5);subplot(2,1,1);plot(yd);ydd=fft(yd,1200);subplot(2,1,2);plot(f(1:600),abs(ydd(1:600))); %sound(yd,FS,bits)Blackman加窗高通滤波器:Fs=22050;[x,FS,bits]=wavread('G:\Users\DP\Desktop\SoundTest.wav'); x=x(:,1);figure(1);subplot(2,1,1);plot(x);%sound(x,FS,bits); %回放语音title('语音信号时域波形图')y=fft(x,3260);f=(FS/1630)*[1:1630];subplot(2,1,2);plot(f(1:1630),abs(y(1:1630)));title('语音信号频谱图');%产生噪声信号并加到语音信号%产生噪声信号并加到语音信号t=0:length(x)-1;zs0=0.05*cos(2*pi*100*t/22050);figure(2);subplot(2,1,1)plot(zs0)title('噪声信号波形');zs1=fft(zs0,1200);%sound(zs,FS,bits); %回放噪音subplot(2,1,2)plot(f(1:600),abs(zs1(1:600)));title('噪声信号频谱');x1=x+zs0';%sound(x1,FS,bits); %回放加入噪声后的语音y1=fft(x1,1200);figure(3);subplot(2,1,1);plot(x1);title('加入噪声后的信号波形'); subplot(2,1,2);plot(f(1:600),abs(y1(1:600)));title('加入噪声后的信号频谱');%高通加窗滤波fp=600,fc=400;wp=2*pi*fp/FS;ws=2*pi*fc/FS;Bt=wp-ws;N0=ceil(11*pi/Bt);N=N0+mod(N0+1,2);wc=(wp+ws)/2/pi;hn=fir1(N-1,wc,'high',blackman(N)); X=conv(hn,x);%sound(X,FS,bits);X1=fft(X,1200);figure(4);subplot(211);plot(X);title('滤波后的信号波形');subplot(212);plot(f(1:600),abs(X1(1:600)));title('滤波后的信号频谱')。
%写上标题%设计低通滤波器:[N,Wc]=buttord()%估算得到Butterworth低通滤波器的最小阶数N和3dB截止频率Wc[a,b]=butter(N,Wc); %设计Butterworth低通滤波器[h,f]=freqz(); %求数字低通滤波器的频率响应figure(2); % 打开窗口2subplot(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('Time waveform');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('butter bandpass filter');grid; %绘制带网格的图像sf=filter(a,b,s); %叠加函数S经过带通滤波器以后的新函数subplot(222);plot(t,sf); %绘制叠加函数S经过带通滤波器以后的时域图形xlabel('Time(seconds)');ylabel('Time waveform');SF=fft(); %对叠加函数S经过带通滤波器以后的新函数进行256点的基—2快速傅立叶变换w=( %新信号角频率subplot(223);plot(')); %绘制叠加函数S经过带通滤波器以后的频谱图title('带通滤波后的频谱图');实例应用:matlab设计的带通滤波器方法改变参数就行了cheb1% wp1=0.45*pi;wp2=0.65*pi;ws1=0.3*pi;ws2=0.75*pi;Rp=1;Rs=40% =============双线型变换法=========================================wp1=0.45*pi; wp2=0.65*pi;ws1=0.3*pi; ws2=0.75*pi;Rp=1; Rs=40;Wp1=tan(wp1/2); Wp2=tan(wp2/2);Ws1=tan(ws1/2); Ws2=tan(ws2/2);BW=Wp2-Wp1; W0=Wp1*Wp2; W00=sqrt(W0);WP=1; WS=WP*(W0^2-Ws1^2)/(Ws1*BW);[N,Wn]=cheb1ord(WP,WS,Rp,Rs,'s');[B,A]=cheby1(N,Rp,Wn,'s');[BT,A T]=lp2bp(B,A,W00,BW);[num,den]=bilinear(BT,A T,0.5);[h,omega]=freqz(num,den,64);subplot(2,2,1);stem(omega/pi,abs(h));xlabel('\omega/\pi');ylabel('|H(z)|');subplot(2,2,2);stem(omega/pi,20*log10(abs(h)));xlabel('\omega/\pi');ylabel('增益.dB');%=============双线性变化法2=================================ws=0.6*pi;Ap=1;As=25;wp=0.4*pi;T=0.001;Fs=1/T;wp=wp/T;ws=ws/T;N=buttord(wp,ws,Ap,As,'s');fprintf('\nN=%d\n',N);wc=wp/((10^(0.1*Ap)-1)^(1/2/N));fprintf('\nwc=%.4e\n',wc);[numa,dena]=butter(N,wc,'s');fprintf('\n');disp('A numerator polynomial');fprintf('%.4e\n',numa);fprintf('\n');disp('A Denominator polynomial');fprintf('%.4e\n',dena);[numd,dend]=bilinear(numa,dena,Fs);w=linspace(0,pi,512);h=freqz(numd,dend,w);norm=max(abs(h));numd1=abs(h)/norm;plot(w/pi,20*log10(numd1));grid;%xlable('Normalized frequency');%ylable('Gain,dB');fprintf('\n');disp('D numerator polynomial');fprintf('%.4e\n',numd);fprintf('\n');disp('D numerator polynomial');fprintf('%.4e\n',dend);% =============直接法================================= wp1=0.45*pi; wp2=0.65*pi;ws1=0.3*pi; ws2=0.75*pi;Rp=1; Rs=40;Wp=[wp1/pi,wp2/pi]; Ws=[ws1/pi,ws2/pi];[N,Wn]=cheb1ord(Wp,Ws,Rp,Rs);[B,A]=cheby1(N,Rp,Wn);[h,omega]=freqz(B,A,64);subplot(2,2,3);stem(omega/pi,abs(h));xlabel('\omega/\pi');ylabel('|H(z)|');subplot(2,2,4);stem(omega/pi,20*log10(abs(h)));xlabel('\omega/\pi');ylabel('增益.dB');%cheby2%% wp1=0.45*pi;wp2=0.65*pi;ws1=0.3*pi;ws2=0.75*pi;Rp=1;Rs=40% =============脉冲响应不变法的程序================================= wp=0.4*pi;ws=0.6*pi;Ap=1;As=25;wp=0.4*pi;T=0.001;Fs=1/T;wp=wp/T;ws=ws/T;N=buttord(wp,ws,Ap,As,'s');fprintf('\nN=%d\n',N);wc=wp/((10^(0.1*Ap)-1)^(1/2/N));fprintf('\nwc=%.4e\n',wc);[numa,dena]=butter(N,wc,'s');fprintf('\n');disp('A numerator polynomial');fprintf('%.4e\n',numa);fprintf('\n');disp('A Denominator polynomial');fprintf('%.4e\n',dena);[numd,dend]=impinvar(numa,dena,Fs);w=linspace(0,pi,512);h=freqz(numd,dend,w);norm=max(abs(h));numd1=abs(h)/norm;plot(w/pi,-20*log10(numd1));grid;%xlable('Normalized frequency');%ylable('Gain,dB');fprintf('\n');disp('D numerator polynomial');fprintf('%.4e\n',numd);fprintf('\n');disp('D numerator polynomial');fprintf('%.4e\n',dend);各种滤波器形式的MATLAB参考程序和仿真内容%*******************************************************************%%mode: 1--巴特沃斯低通;2--巴特沃斯高通;3--巴特沃斯带通;4--巴特沃斯带阻% 5--契比雪夫低通;6--契比雪夫高通;7--契比雪夫带通;8--契比雪夫带阻%fp1,fp2:通带截止频率,当高通或低通时只有fp1有效%fs1, fs2:阻带截止频率,当高通或低通时只有fs1有效%rp: 通带波纹系数%as: 阻带衰减系数%sample: 采样率%h: 返回设计好的滤波器系数%*******************************************************************% function[b,a]=iirfilt(mode,fp1,fp2,fs1,fs2,rp,as,sample)wp1=2*fp1/sample;wp2=2*fp2/sample;ws1=2*fs1/sample;ws2=2*fs2/sample;%得到巴特沃斯滤波器的最小阶数N和3bd频率wnif mode<3[N,wn]=buttord(wp1,ws1,rp,as);elseif mode<5[N,wn]=buttord([wp1 wp2],[ws1 ws2],rp,as);%得到契比雪夫滤波器的最小阶数N和3bd频率wnelseif mode<7[N,wn]=cheb1ord(wp1,ws1,rp,as);else[N,wn]=cheblord([wp1 wp2],[ws1 ws2],rp,as);end%得到滤波器系数的分子b和分母aif mode= =1[b,a]=butter(N,wn);endif mode= =2[b,a]=butter(N,wn,/high/);endif mode= =3[b,a]=butter(N,wn);endif mode= =4[b,a]=butter(N,wn,/stop/);endif mode= =5[b,a]=cheby1(N,rp,wn);endif mode= =6[b,a]=cheby1(N,rp,wn,/high/);endif mode= =7[b,a]=cheby1(N,rp,wn);endif mode= =8[b,a]=cheby1(N,rp,wn,/stop/);endset(gcf,/menubar/,menubar);freq_response=freqz(b,a);magnitude=20*log10(abs(freq_response));m=0:511;f=m*sample/(2*511);subplot(3,1,1);plot(f,magnitude);grid; %幅频特性axis([0 sample/2 1.1*min(magnitude) 1.1*max(magnitude)]); ylabel('Magnitude');xlabel('Frequency-->');phase=angle(freq_response);subplot(3,1,2);plot(f,phase);grid; %相频特性axis([0 sample/2 1.1*min(phase) 1.1*max(phase)]);ylabel('Phase');xlabel('Frequency-->');h=impz(b,a,32); %32点的单位函数响应t=1:32;subplot(3,1,3);stem(t,h);grid;axis([0 32 1.2*min(h) 1.1*max(h)]);ylabel('h(n)');xlabel('n-->');。
各类滤波器的MATLAB程序一、理想低通滤波器IA=imread('lena.bmp');[f1,f2]=freqspace(size(IA),'meshgrid');Hd=ones(size(IA));r=sqrt(f1.^2+f2.^2);Hd(r>0.2)=0;Y=fft2(double(IA));Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA));subplot(2,2,2),imshow(uint8(Ia));figuresurf(Hd,'Facecolor','interp','Edgecolor','none','Facelighting','phong');二、理想高通滤波器IA=imread('lena.bmp');[f1,f2]=freqspace(size(IA),'meshgrid');Hd=ones(size(IA));r=sqrt(f1.^2+f2.^2);Hd(r<0.2)=0;Y=fft2(double(IA));Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya));figuresubplot(2,2,1),imshow(uint8(IA));subplot(2,2,2),imshow(uint8(Ia));figuresurf(Hd,'Facecolor','interp','Edgecolor','none','Facelighting','phong');三、B utterworth低通滤波器IA=imread('lena.bmp');[f1,f2]=freqspace(size(IA),'meshgrid');D=0.3;r=f1.^2+f2.^2;n=4;for i=1:size(IA,1)for j=1:size(IA,2)t=r(i,j)/(D*D);Hd(i,j)=1/(t^n+1);endendY=fft2(double(IA));Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya));figuresubplot(2,2,1),imshow(uint8(IA));subplot(2,2,2),imshow(uint8(Ia));figuresurf(Hd,'Facecolor','interp','Edgecolor','none','Facelighting','phong');四、B utterworth高通滤波器IA=imread('lena.bmp');[f1,f2]=freqspace(size(IA),'meshgrid');D=0.3;r=f1.^2+f2.^2;n=4;for i=1:size(IA,1)for j=1:size(IA,2)t=(D*D)/r(i,j);Hd(i,j)=1/(t^n+1);endendY=fft2(double(IA));Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya));figuresubplot(2,2,1),imshow(uint8(IA));subplot(2,2,2),imshow(uint8(Ia));figuresurf(Hd,'Facecolor','interp','Edgecolor','none','Facelighting','phong');五、高斯低通滤波器IA=imread('lena.bmp');IB=imread('babarra.bmp');[f1,f2]=freqspace(size(IA),'meshgrid');D=100/size(IA,1);r=f1.^2+f2.^2;Hd=ones(size(IA));for i=1:size(IA,1)for j=1:size(IA,2)t=r(i,j)/(D*D);Hd(i,j)=exp(-t);endendY=fft2(double(IA));Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya));figuresubplot(2,2,1),imshow(uint8(IA));subplot(2,2,2),imshow(uint8(Ia));figuresurf(Hd,'Facecolor','interp','Edgecolor','none','Facelighting','phong');六、高斯高通滤波器IA=imread('lena.bmp');IB=imread('babarra.bmp');[f1,f2]=freqspace(size(IA),'meshgrid');%D=100/size(IA,1);D=0.3;r=f1.^2+f2.^2;for i=1:size(IA,1)for j=1:size(IA,2)t=r(i,j)/(D*D);Hd(i,j)=1-exp(-t);endendY=fft2(double(IA));Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya));figuresubplot(2,2,1),imshow(uint8(IA));subplot(2,2,2),imshow(uint8(Ia));figuresurf(Hd,'Facecolor','interp','Edgecolor','none','Facelighting','phong');七、梯形低通滤波器IA=imread('lena.bmp');IB=imread('babarra.bmp');[f1,f2]=freqspace(size(IA),'meshgrid');%D=100/size(IA,1);D0=0.1;D1=0.4;r=sqrt(f1.^2+f2.^2);Hd=zeros(size(IA));Hd(r<D0)=1;for i=1:size(IA,1)for j=1:size(IA,2)if r(i,j)>=D0 & r(i,j)<=D1Hd(i,j)=(D1-r(i,j))/(D1-D0);endendendY=fft2(double(IA));Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya));figuresubplot(2,2,1),imshow(uint8(IA));subplot(2,2,2),imshow(uint8(Ia));figuresurf(Hd,'Facecolor','interp','Edgecolor','none','Facelighting','phong');八、梯形高通滤波器IA=imread('lena.bmp');IB=imread('babarra.bmp');[f1,f2]=freqspace(size(IA),'meshgrid');%D=100/size(IA,1);D0=0.1;D1=0.4;r=sqrt(f1.^2+f2.^2);Hd=ones(size(IA));Hd(r<D1)=0;for i=1:size(IA,1)for j=1:size(IA,2)if r(i,j)>=D0 & r(i,j)<=D1Hd(i,j)=(D0-r(i,j))/(D0-D1);endendendY=fft2(double(IA));Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya));figuresubplot(2,2,1),imshow(uint8(IA));subplot(2,2,2),imshow(uint8(Ia));figuresurf(Hd,'Facecolor','interp','Edgecolor','none','Facelighting','phong');九、用其他方法编写的理想低通、理想高通、Butterworth低通、同态滤波程序1、理想低通i1=imread('lena.bmp');i2=imnoise(i1,'salt & pepper',0.1);f=double(i2);k=fft2(f);g=fftshift(k);[N1,N2]=size(g);d0=50;u0=floor(N1/2)+1;v0=floor(N2/2)+1;for i=1:N1for j=1:N2d=sqrt((i-u0)^2+(j-v0)^2);if d<=d0h=1;elseh=0;endy(i,j)=g(i,j)*h;endendy=ifftshift(y);E1=ifft2(y);E2=real(E1);figuresubplot(2,2,1),imshow(uint8(i1));subplot(2,2,2),imshow(uint8(i2));subplot(2,2,3),imshow(uint8(E2));2、理想高通i1=imread('lena.bmp');i2=imnoise(i1,'salt & pepper',0.1); f=double(i2);k=fft2(f);g=fftshift(k);[N1,N2]=size(g);n=2;d0=10;u0=floor(N1/2)+1;v0=floor(N2/2)+1;for i=1:N1for j=1:N2d=sqrt((i-u0)^2+(j-v0)^2);if d<=d0h=0;else h=1;endy(i,j)=g(i,j)*h;endendy=ifftshift(y);E1=ifft2(y);E2=real(E1);figuresubplot(2,2,1),imshow(uint8(i1)); subplot(2,2,2),imshow(uint8(i2)); subplot(2,2,3),imshow(uint8(E2));3、Butterworth低通i1=imread('lena.bmp');i2=imnoise(i1,'salt & pepper',0.1); f=double(i2);k=fft2(f);g=fftshift(k);[N1,N2]=size(g);n=2;d0=50;u0=floor(N1/2)+1;v0=floor(N2/2)+1;for i=1:N1for j=1:N2d=sqrt((i-u0)^2+(j-v0)^2);h=1/(1+(d/d0)^(2*n));y(i,j)=g(i,j)*h;endendy=ifftshift(y);E1=ifft2(y);E2=real(E1);figuresubplot(2,2,1),imshow(uint8(i1)); subplot(2,2,2),imshow(uint8(i2)); subplot(2,2,3),imshow(uint8(E2));4、同态滤波I=rgb2gray(imread('fabric00.bmp')); [M,N]=size(I);T=double(I);L=log(T);F=fft2(L);A=2;B=0.3;for i=1:Mfor j=1:ND(i,j)=((i-M/2)^2+(j-N/2)^2);endendc=1.1;%锐化参数D0=max(M,N);H=(A-B)*(1-exp(c*(-D/(D0^2))))+B;F=F.*H;F=ifft2(F);Y=exp(F);figuresubplot(1,2,1),imshow(I);subplot(1,2,2),imshow(uint8(real(Y)));十、G abor滤波器。