应用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);

end

end %将M*N维矩阵变成1维矩阵

fs=1000;%扫描频率1kHz

T=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);

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

x7(i,j)=x6(N*(i-1)+j);

end

end

%%设计FIR滤波器并分析相关指标

wc=280*2/fs; %6dB截止频率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('原始灰度图像幅度谱');