应用Matlab对图像信号进行频谱分析及滤波
应用Matlab对图像信号进行频谱分析及滤波实验目的1. 巩固所学的数字信号处理理论知识，理解信号的采集、处理、传输、显示和存储过程；2. 综合运用专业及基础知识，解决实际工程技术问题的能力；3. 学习资料的收集与整理，学会撰写课程设计报告。实验环境1. 微型电子计算机（PC）；2. 安装Windows 10操作系统，MATLAB7.0，Formatfactory，绘画板等开发工具。实验原理 在Matlab软件平台下，读取和显示彩色图像数据的相关函数和调用方法如下： 实验内容和任务要求1. 选取一张彩色图像（注意不能出现雷同，否则记为0分），提取图像的灰度值，并显示出灰度图像。2. 在图像中增加正弦噪声信号(自己设置几个频率的正弦信号)，显示出加入噪声信号后的灰度图像。3. 给定滤波器的性能指标，分别设计FIR和IIR数字滤波器，并画出滤波器的幅频响应曲线。4. 用自己设计的滤波器对含噪声图像信号进行滤波，显示出滤波后的灰度图像。5. 对原始灰度图像、加入噪声信号的灰度图像和滤波后的灰度图像进行频谱分析和对比，分析信号的变化。实验分析本实验要求用Matlab软件完成对图像信息的处理。对于任务1，这里采用了一张jpg格式的张学友新专辑醒着做梦的封面图片，保存在Matlab的work文件夹下。采用imread()函数读取，并利用rgb2gray()函数将其转化为二维的灰度图像（原始的数据类型是unit8型，需要将其转化为可用于计算的double类型），并利用imshow()函数将其显示出来；对于任务2，在加入噪声前，需要先将二维数据利用循环嵌套语句转化为一维数据，然后加入三个高频噪声，再利用循环嵌套语句转化为二维的数据，利用imshow()函数显示出来；对于任务3，这里分别设计了满足一定指标的IIR低通滤波器（巴特沃斯）和FIR低通滤波器（哈明窗）并对其相关指标进行了分析。对于任务4，利用任务3中设计好的两个滤波器分别对加噪后的灰度图像进行滤波（filter()函数），并分别显示滤波后的灰度图像；对于任务5，利用快速傅里叶变换算法（FFT）对各阶段数据分别进行频谱分析，并将它们的频谱绘制在同一张图上作为对比。Matlab代码clc;close all;clear all;%图像的读取以及转换x=imread('hh.jpg');%读取jpg图像x1=rgb2gray(x);%生成M*N的灰度图像矩阵M,N=size(x1);%求图像规模%生成原始序列并求频率响应x2=im2double(x1);x3=zeros(1,M*N);%初始化for i=1:M for j=1:N x3(N*(i-1)+j)=x2(i,j); endend %将M*N维矩阵变成1维矩阵fs=1000;%扫描频率1kHzT=1/fs; %扫描时间间隔L=length(x3);%计算序列长度n=0:L-1;Xk3=fft(x3);%快速傅里叶变换Xf3=fftshift(Xk3);%中心对称变换f=(n/L-1/2)*fs;%横坐标变幻%生成带有噪声的序列并求频率响应fz1=356;fz2=383;fz3=427;%三个噪声频率xz=0.4*sin(2*pi*fz1*n*T)+0.7*sin(2*pi*fz2*n*T)+0.5*sin(2*pi*fz3*n*T);%噪声序列x4=x3+xz;%加入噪声的序列x5=zeros(M,N);for i=1:M for j=1:N x5(i,j)=x4(N*(i-1)+j); endend %一维变M*N矩阵Xk4=fft(x4);Xf4=fftshift(Xk4);%设计IIR滤波器并分析相关指标wp=250*2/fs;ws=300*2/fs;Rp=3;Rs=20;%通带截止频率250Hz,阻带截止频率300Hz,通带衰减3dB,阻带衰减20dBNm,Wc=buttord(wp,ws,Rp,Rs);%计算满足指标的最小阶数以及3dB截止频率b,a=butter(Nm,Wc); %计算滤波器的分子分母系数H=freqz(b,a,f*2*pi/fs);%计算滤波器频率响应mag=abs(H);pha=angle(H);%幅度响应和相位响应mag1=20*log(mag+eps)/max(mag); %将幅频响应转化为dB的形式%用IIR滤波器对带噪序列进行滤波并求频率响应x6=filter(b,a,x4);%用IIR滤波Xk6=fft(x6);Xf6=fftshift(Xk6);x7=zeros(M,N);for i=1:M for j=1:N x7(i,j)=x6(N*(i-1)+j); endend%设计FIR滤波器并分析相关指标wc=280*2/fs; %6dB截止频率280kHzfx=0 wc wc 1;m=1 1 0 0;%理想频幅响应b1=fir2(40,fx,m,hamming(41);%计算FIR滤波器多项式系数（不声明窗默认为Hamming窗）H1=freqz(b1,1,f*2*pi/fs);%计算滤波器频率响应mag2=abs(H1);pha1=angle(H1);%幅度响应和相位响应mag3=20*log(mag2+eps)/max(mag2); %将幅频响应转化为dB的形式%用FIR滤波器对带噪序列进行滤波并求频率响应x8=filter(b1,1,x4);%用FIR进行滤波grd=grpdelay(b1,1,f*2*pi/fs);%计算群延时K=round(grd(1);%修正因子(冒号操作做索引时，需要整型数操作)x8=x8(K+1):L),x8(1:K);%对群延迟进行修正Xk8=fft(x8);Xf8=fftshift(Xk8);x9=zeros(M,N);for i=1:M for j=1:N x9(i,j)=x8(N*(i-1)+j); endend%绘制图像figure(1);subplot(2,2,1);imshow(x2);title('原始灰度图像');subplot(2,2,2);imshow(x5);title('加入噪声后灰度图像');subplot(2,2,3);imshow(x7);title('IIR滤波器滤波后灰度图像');subplot(2,2,4);imshow(x9);title('FIR滤波器滤波后灰度图像');figure(2);subplot(4,1,1);plot(f,abs(Xf3)*2/L,'r-');title('原始灰度图像幅度谱');subplot(4,1,2),plot(f,abs(Xf4)*2/L,'r-');title('加入噪声后灰度图像幅度谱');subplot(4,1,3);plot(f,abs(Xf6)*2/L,'r-');title('IIR滤波器滤波后灰度图像幅度谱');subplot(4,1,4);plot(f,abs(Xf8)*2/L,'r-');title('FIR滤波器滤波后灰度图像幅度谱');figure(3);subplot(3,2,1);plot(f,mag);grid;title('IIR滤波器幅度响应');xlabel('f/Hz');ylabel('幅度');subplot(3,2,2);plot(f,mag2);grid;title('FIR滤波器幅度响应');xlabel('f/Hz');ylabel('幅度');subplot(3,2,3);plot(f,pha);grid;title('IIR滤波器相位响应');xlabel('f/Hz');ylabel('相位');subplot(3,2,4);plot(f,pha1);grid;title('FIR滤波器相位响应');xlabel('f/Hz');ylabel('相位');subplot(3,2,5);plot(f,mag1);grid;title('IIR滤波器幅度响应(dB)');xlabel('f/Hz');ylabel('幅度/dB');subplot(3,2,6);plot(f,mag3);grid;title('FIR滤波器幅度响应(dB)');xlabel('f/Hz');ylabel('幅度/dB');实验结果及分析从实验结果来看，滤波效果还是相当不错的。在满足相同的指标情况下FIR滤波器所需要的阶数远高于IIR滤波器，但是FIR滤波器的相位响应是线性的，滤波器通带群延时是常数，而IIR滤波器是非线性的，滤波器通带群延时非常数。FIR滤波器产生的群延迟样本数可以计算出来，并进行调整。IIR滤波器则不可以。调整代码为：grd=grpdelay(b1,1,f*2*pi/fs);%计算群延时K=round(grd(1);%修正因子(冒号操作做索引时，需要整型数操作)x8=x8(K+1):L),x8(1:K);%对群延迟进行修正FIR滤波器滤波后图像群延迟调整前后对比：将两种滤波器滤波后的图片放大后对比：对比可发现IIR滤波器滤波后使图片产生了重影，而FIR滤波器没有明显失真。对于语音系统，対相位要求不是主要的，因此，选用IIR滤波器较为合适，可以充分发挥其经济和高效的特点；图像信号和数据传输等以波形携带信息的系统对线性相位要求较高，因此采用FIR滤波器较好。实验总结通过本次实验，加深了我们对信号频谱分析和数字滤波器设计的知识的理解，并对滤波器有了更进一步的认识，掌握了利用滤波器处理图片的方法，理解了设计指标的工程概念，认识到了不同类型滤波器的特性和适用范围。实验中，我们对Matlab中一些函数的用法还不清楚，后面经过查资料以及尝试后均得到了解决。总的来说，本次实验基本上达到了预期的实验效果，是一次比较成功的实验。8
