写在前面,这是一篇有关于"XX"大学光学"X"程的数值"XX"课作业的讲解汇总,江月年年望相似,作业也是如此。
本人中秋国庆回不了家,无事可干,故于2023年9月30日决定写下该回答,仅供参考,节约诸位时间。如果Anybody因为照抄导致你自己挂科,并且让老师抗拒这种形式的交流,那你是真该死
最后,愿本回答会成为帮助各位师弟、师妹的渡河之桥,更快抵达你的诗和远方
第一部分:实现的功能
A.要求功能一
- 读入一幅灰度级不协调的图像,如pout.tif
a.分析其直方图(文字说明);
b.设计灰度变换公式;
c.调用imadjust函数,调整变换表达式参数到图像的灰度级分布均衡。 - 对自选一幅图像调价噪声,使用均值滤波和中值滤波去噪,并显示和保存
B.要求功能二
- 用傅立叶变换求由两个时间周期分别为0.2 和5的正弦函数;函数,y1,y2之和组成的函数的频谐。y1和y2的振幅分别为5和1。
y1=5 * sin(2 * pi * 5 * x)-----y2 = sin(2 * pi * 0.2 * x)------y= y1+y2 - a.在rice.tif中加噪声,并用2阶Butterworth低通滤波器去除噪声.
b.用高斯高通滤波器对rice.tif原图图像滤波
第二部分:基本说明
A.用到的函数
imhist()
基本用法:用于显示直方图,要创建显示窗口
figure('Name','直方图'); %创建显示直方图的窗口
imhist(pic); %显示直方图
fspecial()
基本用法:fspecial(‘要的滤波器类型’,参数)
fil = fspecial('average',[3,3]); %创建3 * 3的均值滤波矩阵。
滤波器字符 | 滤波器类型 |
---|---|
‘average’ | 平均值滤波器 |
‘disk’ | 圆形平均值滤波器 (pillbox) |
‘gaussian’ | 高斯低通滤波器。不推荐。请改用 imgaussfilt 或 imgaussfilt3。 |
‘laplacian’ | 逼近二维拉普拉斯算子 |
‘log’ | 高斯拉普拉斯滤波器 |
‘motion’ | 逼近相机的线性运动 |
‘prewitt’ | Prewitt 水平边缘强调滤波器 |
‘sobel’ | Sobel 水平边缘强调滤波器 |
imfilter()
基本用法:imfilter函数叫做实现线性空间滤波函数,功能是对任意类型数组或多维图像进行滤波,函数形式是 pic = imfilter(原图片矩阵,滤波矩阵)。
fil = fspecial('average',[3,3]); %创建3 * 3的滤波矩阵。
pic_ave = imfilter(pic_noise,fil,'symmetric'); %均值滤波,边界之外的输入数组值是通过沿数组边界对数组进行镜面反射得到
imadjust()
基本用法:imadjust(f,[low_in high_in],[low_out high_out],gamma)
pic_ad = imadjust(pic,[0.29,0.63],[0,255]);
%调整矩阵灰度范围,获得输入图片0.29*255到0.63*255的灰度值,线性扩展变为灰度值0到255的输出图
fft()、fft2()
基本用法:傅里叶变换函数,就是简单的直接用
y = [1,2,3,4,5,6,7,8,9,10]; %原函数,或者说一阶的行列式
Fy = fft(y); %对原函数进行1阶快速傅里叶变换
pic = [1,2,3;3,2,1;6,6,6]; %2阶矩阵
FFT_pic = fft2(pic); %对原图矩阵进行2阶离散傅里叶变换
结果展示:
% y = [1,2,3,4,5,6,7,8,9,10];
Fy =
21.0000 + 0.0000i -3.0000 + 5.1962i -3.0000 + 1.7321i -3.0000 + 0.0000i -3.0000 - 1.7321i -3.0000 - 5.1962i
% pic = [1,2,3;3,2,1;6,6,6];
FFT_pic =
30.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i
-6.0000 +10.3923i -3.0000 + 0.0000i -1.5000 - 2.5981i
-6.0000 -10.3923i -1.5000 + 2.5981i -3.0000 - 0.0000i
ifft()、ifft2()
基本用法:傅里叶逆变换函数,直接用就好了,同上
fftshift()、iffshift()
基本用法:它本质就是一个矩阵交换的方法,目的是将零频量移到中间。举个例子:
%ffftshift和ifftshift的例子
% 元素个数为偶数
x=[1,2,3,4,5,6];
x_fftshift=fftshift(x)
x_ifftshift=ifftshift(x)
% 元素个数为奇数
y=[1,2,3,4,5];
y_fftshift=fftshift(y)
y_ifftshift=ifftshift(y)
运行结果
%运行结果为
%######## x=[1,2,3,4,5,6]; 元素个数为偶数
x_fftshift =
4 5 6 1 2 3
x_ifftshift =
4 5 6 1 2 3
%######## y=[1,2,3,4,5]; 元素个数为奇数
y_fftshift =
4 5 1 2 3
y_ifftshift =
3 4 5 1 2
B.难点思路分析
首先,我觉得这东西没啥难点,和同学交流,发现他们对时间周期,采样间隔,采样频率之类的概念有巨多的疑问,所以我就先说说这个问题
1、概念理解
首先,你要知道对于计算机,没有常规上的变量而言,有的只是一个个有对应关系的值。对于一个sin(2 * pi * x)函数的例子,获得的永远是一个个有一定对应关系的值,每一个 t 的值确定下来,y 就确定了
x = 0:0.1:1
fx = sin(2*pi*x)
figure;
plot(x,y,'r*');
% 运行结果展示:
x =
0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000
fx =
0 0.5878 0.9511 0.9511 0.5878 0.0000 -0.5878 -0.9511 -0.9511 -0.5878 -0.0000
也就是说,傅里叶变换中的 f(x) 其实是一个个数值,x 是从 0 到 N-1的一个个整数。带入一个 u 的值,将每一个 x 和 f(x) 累加起来,就会得到一个 F(u) 的值。
计算机中的傅里叶变换,要注意的是,f(x) 中的 x 对matlab来说,其实是数组的第几个元素。 也就是说 x 一定是从 1 到 N 。与你取的 x 值毫无关系
这是计算机进行傅里叶变换的基本逻辑,所以怎么取好 x 的值很重要,主要表现在对 fx 数组值的影响上。
A.采样间隔(核心)
对于上面sin(2 * pi * t)函数的例子,我们获得 ft 数组,对其傅里叶变换最重要的是, ft 数组中的值有没有涵盖sin函数中有特征的值,ft 数值的值。下图展示一个典型的采样间隔错误
对于sin(2 * pi * t)我们 t 的采样间隔至少要大于0.5,这样才能大概取得函数中有鲜明特征的值
采样频率 = 1 / 采样间隔;
采样数 = 周期长度 / 采样间隔 + 1;(加1 是加上开始的原点,即 100个间隔,要采101个点)
2、滤波器去噪声的实现
1、傅里叶变换,获得频谱
2、用滤波器作用频谱,即滤波器乘以频谱,获得滤波后频谱
3、对滤波后频谱进行傅里叶逆变换
A.循环
我交的作业就是用这个,因为这个会提高代码的可读性,你去看下我的作业第二题,你就知道了,就不重复思路代码了。
B.矩阵基本积(哈达玛积)
% 在rice.tif中加噪声,并用2阶Butterworth低通滤波器去除噪声.
pic_path = 'rice.png' ; %输入图片路径
pic = imread(pic_path); %打开图片
pic=im2double(pic); %图片类转为double类
figure(); %创建显示窗口
subplot(1,2,1); %组合绘图1
title('原图'); %原图标题
imshow(pic); %显示图片
%######### 在rice.tif中加噪声 ###############
pic_noise = imnoise(pic,'salt & pepper'); %添加椒盐噪声。高斯噪声'gaussian';泊松噪声'poisson';高斯白噪声'localvar';乘性噪声'speckle'。
subplot(1,2,2); %组合绘图2
title('噪声图'); %噪声图标题
imshow(pic_noise); %显示图片
%######### 傅里叶变换 ############
FFT_pic = fft2(pic); %原图矩阵进行2阶傅里叶变换
FFT_pic_shift = fftshift(fft2(pic)); %变换后矩阵进行频移
[u_max,v_max] = size(FFT_pic_shift); %获取变换后矩阵的行列数
%######### 定义频谱中的参量u、v和中心频率 ##########
n = 2; %Butterworth的阶数
D0 = 20; %Butterworth的截止频率
u0 = u_max / 2; %获取中心频率
v0 = v_max / 2; %获取中心频率
u = 1 : u_max; %定义 FFT_pic(u,v) 中的频率 u 的行列式
v = 1 : v_max; %定义 FFT_pic(u,v) 中的频率 v 的行列式
%######## 矩阵基本积思路 ##################
D=sqrt( (u - u0) .^ 2 + (v - v0) .^ 2); %定义一个频域滤波的范围
H = 1 ./ (1 + (D ./ D0) .^ (2*n)); %定义一个对 频域滤波范围内外 的作用方法--即--滤波器
FFT_pic_fre = fftshift(FFT_pic_shift .* H); %用滤波器作用 频移后的频谱 ,并再次进行频移
pic_butter = abs(ifft2(FFT_pic_fre)); %对矩阵进行逆变换后求绝对值
figure('Name','变换后butter降噪图'); %创建显示变换后降噪图的窗口
imshow(pic_butter); %显示图片
第三部分:我的作业
第一题
% 文件名: DSH_2_1.m
% 版本: 1.0
% 创建时间: 2023年09月19日 星期二 18:37
% 功能: 1. 读入一幅灰度级不协调的图像,如pout.tif
% a.分析其直方图(文字说明);b.设计灰度变换公式;c.调用imadjust函数,调整变换表达式参数到图像的灰度级分布均衡。
% 2.对自选一幅图像调价噪声,使用均值滤波和中值滤波去噪,并显示和保存
%功能: 1. 读入一幅灰度级不协调的图像,如pout.tif
pic_path = 'pout.tif' ; %输入图片路径
pic = imread(pic_path); %打开图片
figure('Name','原图'); %创建显示原图的窗口
imshow(pic); %显示图片
figure('Name','直方图'); %创建显示直方图的窗口
imhist(pic); %显示直方图
%分析:直方图横坐标为灰度,纵坐标为灰度点个数,灰度大多数集中在75-175,所以我们要让其灰度范围变广,不那么集中。此外图片的灰度变化不明显,人眼不能很好的区分,要让各灰度间隔变大。
%设计公式:公式为线性变化,灰度范围由[75,160]线性变化到[0,255],变换公式应该为 y = (255 / (160 - 75)) * (x - 75),化简为y = 3 * x -255 。
pic_ad = imadjust(pic,[0.29,0.63],[]); %调整矩阵灰度范围
figure('Name','调整图'); %创建显示直方图的窗口
imshow(pic_ad); %显示直方图
figure('Name','调整直后方图'); %创建显示调整后直方图的窗口
imhist(pic_ad); %显示直方图
%功能: 2.对自选一幅图像调价噪声,使用均值滤波和中值滤波去噪,并显示和保存
pic_path = 'pout.tif' ; %输入图片路径
win_path = 'C:\Users\DSH\Desktop\'; %输入保存路径--默认是DSH的桌面
pic = imread(pic_path); %打开图片
pic_noise = imnoise(pic,'salt & pepper'); %添加椒盐噪声。高斯噪声'gaussian';泊松噪声'poisson';高斯白噪声'localvar';乘性噪声'speckle'。
figure('Name','噪声图'); %创建显示噪声图的窗口
imshow(pic_noise); %显示图片
fil = fspecial('average',[3,3]); %创建3 * 3的滤波矩阵。
pic_ave = imfilter(pic_noise,fil,'symmetric'); %均值滤波,边界之外的输入数组值是通过沿数组边界对数组进行镜面反射得到
figure('Name','均值降噪图'); %创建显示均值降噪图的窗口
imshow(pic_ave); %显示图片
imwrite(pic_ave,strcat(win_path,'pic_ave.tif')); %保存均值降噪图到桌面
pic_med = medfilt2(pic_noise,[3,3]); %中值滤波,3 * 3范围选一个
figure('Name','中值降噪图'); %创建显示中值降噪图的窗口
imshow(pic_med); %显示图片
imwrite(pic_med,strcat(win_path,'pic_med.tif')); %保存中值降噪图到桌面
fprintf('您好,世界\t来自DSH的问候\n');
第二题
% 文件名: DSH_2_2.m
% 版本: 1.0
% 创建时间: 2023年09月20日 星期二 16:28
% 功能: 1. 用傅立叶变换求由两个时间周期分别为0.2 和5的正弦函数;函数,y1,y2之和组成的函数的频谐。y1和y2的振幅分别为5和1。
% y1=5 * sin(2 * pi * 5 * x) y2 = sin(2 * pi * 0.2 * x) y= y1+y2
% 2. a.在rice.tif中加噪声,并用2阶Butterworth低通滤波器去除噪声.
% b.用高斯高通滤波器对rice.tif原图图像滤波
% %1、用傅立叶变换求由两个时间周期分别为0.2 和5的正弦函数;函数,y1,y2之和组成的函数的频谐。y1和y2的振幅分别为5和1。 y1=5 * sin(2 * pi * 5 * x) y2 = sin(2 * pi * 0.2 * x) y= y1+y2
Fs = 20; %两函数频率分别为 5 和 0.2 ,采样频率应该大于 2 * 5 = 10 Hz; 这里取 20Hz
dx = 1 / Fs; %采样间隔为采样频率的倒数,即1/20;
L = 201; %采样数乘以采样间隔要大于函数一个周期(分别为 0.2 和 5),即采样数(L)>最大周期(5)/采样间隔(1/20);即 L > 100 ,这里取201
df = Fs / (L - 1); %频率间隔,即画出的频谱图的间隔
x = (0 : (L-1)) * dx; %x的取值从0*dx到(L-1)*dx,即 0,1/20,2/20,3/20,
y1 = 5 * sin(2 * pi * 5 * x); %函数y1
y2 = sin(2 * pi * 0.2 * x); %函数y2
y = y1 + y2; %两函数之和
fy =abs(fftshift(fft(y))); %fft傅里叶变换;之后fftshift将零频分量移动到数组中心,重新排列傅里叶变换;最后取其绝对值。
fx = - Fs / 2 : df :Fs/2; %频谱从-Fs/2到Fs/2,间隔为df
figure('name','频谱图'); %创建显示频谱图的窗口
plot(fx,fy); %显示频谱图
%2. a.在rice.tif中加噪声,并用2阶Butterworth低通滤波器去除噪声. b.用高斯高通滤波器对rice.tif原图图像滤波
pic_path = 'rice.png' ; %输入图片路径
pic = imread(pic_path); %打开图片
pic=im2double(pic); %
figure(); %创建显示窗口
subplot(1,2,1); %组合绘图1
title('原图'); %原图标题
imshow(pic); %显示图片
pic_noise = imnoise(pic,'salt & pepper'); %添加椒盐噪声。高斯噪声'gaussian';泊松噪声'poisson';高斯白噪声'localvar';乘性噪声'speckle'。
subplot(1,2,2); %组合绘图2
title('噪声图'); %噪声图标题
imshow(pic_noise); %显示图片
FFT_pic = fft2(pic); %原图矩阵进行2阶傅里叶变换
FFT_pic_shift = fftshift(fft2(pic)); %变换后矩阵进行频移
[u_max,v_max] = size(FFT_pic_shift); %获取变换后矩阵的行列数
%#########a.在rice.tif中加噪声,并用2阶Butterworth低通滤波器去除噪声######################
n = 2; %Butterworth的阶数
D0 = 20; %Butterworth的截止频率
u0 = u_max / 2; %获取中心频率
v0 = v_max / 2; %获取中心频率
FFT_pic_butter_shift = ones(u_max,v_max); %创建全一矩阵
for x=1:u_max
for y=1:v_max
D=sqrt((x-u0)^2+(y-v0)^2); %计算点(x,y)到中心点的距离
h=1/(1+(D/D0)^(2*n)); %计算巴特沃斯滤波器
FFT_pic_butter_shift(x,y)=FFT_pic_shift(x,y)*h; %用滤波器乘以主函数
end
end
FFT_pic_butter = fftshift(FFT_pic_butter_shift); %对矩阵进行频移
pic_butter = abs(ifft2(FFT_pic_butter)); %对矩阵进行逆变换后求绝对值
figure('Name','变换后butter降噪图'); %创建显示变换后降噪图的窗口
imshow(pic_butter); %显示图片
%#########b.用高斯高通滤波器对rice.tif原图图像滤波#################################
u0 = u_max / 2; %获取中心频率
v0 = v_max / 2; %获取中心频率
D0 = 20; %gauss的截止频率
FFT_pic_gauss_shift = ones(u_max,v_max); %创建全一矩阵
for x=1:u_max
for y=1:v_max
D=sqrt((x-u0)^2+(y-v0)^2); %计算点(x,y)到中心点的距离
h = 1 - exp(-D^2 / (2 * D0 ^ 2)); %计算高斯滤波器
FFT_pic_gauss_shift(x,y)=FFT_pic_shift(x,y)*h; %用滤波器乘以主函数
end
end
FFT_pic_gauss = fftshift(FFT_pic_gauss_shift); %对矩阵进行频移
pic_gauss = abs(ifft2(FFT_pic_gauss)); %对矩阵进行逆变换后求绝对值
figure('Name','变换后gauss降噪图'); %创建显示变换后降噪图的窗口
imshow(pic_gauss); %显示图片
fprintf('您好,世界\t来自DSH的问候\n');