发布时间 : 星期四 文章应用Matlab对图像信号进行频谱分析及滤波更新完毕开始阅读526749c027d3240c8447efbd
应用Matlab对图像信号进行频谱分析及滤波
实验目的
1. 巩固所学的数字信号处理理论知识，理解信号的采集、处理、传输、显示和存
储过程；
2. 综合运用专业及基础知识，解决实际工程技术问题的能力； 3. 学习资料的收集与整理，学会撰写课程设计报告。
实验环境
1. 微型电子计算机（PC）； 2. 安装Windows 10操作系统，MATLAB7.0，Formatfactory，绘画板等开发工具。
实验原理
在Matlab软件平台下，读取和显示彩色图像数据的相关函数和调用方法如下：
实验内容和任务要求
1. 选取一张彩色图像（注意不能出现雷同，否则记为0分），提取图像的灰度值，
并显示出灰度图像。 2. 在图像中增加正弦噪声信号(自己设置几个频率的正弦信号)，显示出加入噪声信
号后的灰度图像。
3. 给定滤波器的性能指标，分别设计FIR和IIR数字滤波器，并画出滤波器的幅
频响应曲线。
4. 用自己设计的滤波器对含噪声图像信号进行滤波，显示出滤波后的灰度图像。 5. 对原始灰度图像、加入噪声信号的灰度图像和滤波后的灰度图像进行频谱分析
和对比，分析信号的变化。
1
实验分析
本实验要求用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); end
end %将M*N维矩阵变成1维矩阵 fs=1000;%扫描频率1kHz T=1/fs; %扫描时间间隔
2
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); end
end %一维变M*N矩阵 Xk4=fft(x4);
Xf4=fftshift(Xk4);
%%设计IIR滤波器并分析相关指标
wp=250*2/fs;ws=300*2/fs;Rp=3;Rs=20;%通带截止频率250Hz,阻带截止频率300Hz,通带衰减3dB,阻带衰减20dB
[Nm,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
3
x7(i,j)=x6(N*(i-1)+j); end end
%%设计FIR滤波器并分析相关指标 wc=280*2/fs; mB截止频率280kHz
fx=[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); end end
%%绘制图像 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('原始灰度图像幅度谱');
4