【信号与系统 x MATLAB】实验三

实验摘要：

用 MATLAB 进行信号与系统课程第三次实验

主要内容有：

图片的高频信息与低频信息（合成、分解图片）
频域制作数字盲水印和去除数字盲水印
傅里叶变换、傅里叶逆变换与其时移特性

实验题目

图片的高频信息与低频信息合成图片。

合成图片。找两张轮廓比较像的图片 A 和 B，有一张是你本人。提取一张照片的低频信息，另一张图片的高频信息，结合这两个照片。设置不同的频率门限，组合照片，组合的效果是，放大看是 A，缩小看是 B。例如以下两张图片

分解图片。把下图爱因斯坦和玛丽莲梦露分开（此图缩小是玛丽莲梦露，截取不同的频段）

频域制作数字盲水印和去除数字盲水印 https://www.zhihu.com/question/50735753，看懂，想想，有想法写出来，做一个好玩的东西。
傅里叶变换两题

使用 2 维离散余弦变换（DCT）分解合成图像的高频和低频

clear

合成图像

RGB1 = imread("doge_young.png");
I1 = im2gray(RGB1);
RGB2 = imread("doge_old.png");
I2 = im2gray(RGB2);
% 使用 dct2 对灰度图像执行 2-D DCT
J1 = dct2(I1);
J2 = dct2(I2);
% 高低频分界线
d = 260;
% 去除低频
J_high = J1;
J_high(abs(J1) > d) = 0;
% 去除高频
J_low = J2;
J_low(abs(J2) <= d) = 0;
% 使用逆 DCT 函数 idct2 重建图像
K = idct2(J_low+J_high);
K = rescale(K);
% 显示图像
montage({I1,I2})
title('原始图像');
montage({K})
title('合并图像');
imwrite(K,"such.png")

分解图像

RGB = imread('m_a.png');
I = im2gray(RGB);
% 使用 dct2 对灰度图像执行 2-D DCT
J = dct2(I);
% 使用对数标度显示变换后的图像
% 注意到大部分能量位于左上角
imshow(log(abs(J)),[])
colormap parula
colorbar
% 高低频分界线
d = 340;
% 去除高频
J_low = J;
J_low(abs(J) <= d) = 0;
% 去除低频
J_high = J;
J_high(abs(J) > d) = 0;
% 使用逆 DCT 函数 idct2 重建图像
K_low = idct2(J_low);
K_low = rescale(K_low);
K_high = idct2(J_high);
K_high = rescale(K_high);
% 并排显示图像
montage({K_low,K_high})
title('低频图像/高频图像');

参考：MATLAB 帮助文档例程：Remove High Frequencies in Image using 2-D DCT

输出

给上一题的图片加上盲水印

以下的MATLAB实现的一个极其简单的盲水印demo：

clear

% 加水印
img = double(imread("such.png"));
wm = double(imread("blind.png"));
imshow(wm);
title('水印图像');
d = 5; % 水印深度
fimgre = fft2(img) + wm*d;
imgre = uint8(real(ifft2(fimgre)));
% 并排显示加水印前后图像
montage({uint8(img),imgre})
title('原始图像/加入水印后的图像');

% 解水印
wm = fft2(imgre) - fft2(img);
imshow(mat2gray(real(wm)/d));
title('解出水印');

输出

傅里叶变换

因为涉及到符号计算，于是直接上 Mathematica 写了：

Mathematica 代码：

$\text{FourierTransform}\left[e^{-2 | t| },t,\omega ,\text{FourierParameters}\to {1,-1}\right]$

$\text{InverseFourierTransform}\left[\frac{1}{\omega ^2+1},\omega ,t,\text{FourierParameters}\to {1,-1}\right]$

$\text{FourierTransform}\left[\frac{1}{2} e^{-2 t} \theta (t),t,\omega ,\text{FourierParameters}\to {1,-1}\right]$

$Plot[{Norm[-(I/(2 (-2 I + [Omega])))],
Arg[-(I/(2 (-2 I + [Omega])))]}, {[Omega], -4.2, 4.2},
PlotRange -> {-4, 4}]$

$FourierTransform[E^(-2*(t - 1))/2*HeavisideTheta[t - 1], t, [Omega],
FourierParameters -> {1, -1}]$

$\text{Plot}\left[\left{\left| -\frac{i \cos (\omega )+\sin (\omega )}{2 (-2 i+\omega )}\right| ,\arg \left(-\frac{\sin (\omega )+i \cos (\omega )}{2 (\omega -2 i)}\right)\right},{\omega ,-4.2,4.2},\text{PlotRange}\to {-4,4}\right]$

输出：

后来有同学问我 MATLAB 怎么写，于是用 MATLAB 又写了一遍~~（MATLAB 符号计算居然也能用，但依然大大不如 Mathematica）~~

clear
syms t w

f = exp(-2*abs(t));
fourier(f)

F = 1/(1+w^2);
ifourier(F)

T = -5:0.1:5;
f1 = exp(-2*t)/2*heaviside(t);
fFwl = matlabFunction(fourier(f1))

plot(T,abs(fFwl(T)),'.')
hold on 
plot(T,angle(fFwl(T)),'.')

% 时域移动t-1
f1 = exp(-2*(t-1))/2*heaviside(t-1);
fFwl = matlabFunction(fourier(f1))

plot(T,abs(fFwl(T)))
plot(T,angle(fFwl(T)))
hold off

grid on