close all
clear
%% 1

%The signal in MATLAB is a discrete signal, and the sinusoidal signal is a continuous signal, so it needs to be sampled to discretize it.
%MATLAB中的信号为离散信号，而正弦信号为 连续信号，所以需要 采样 将其 离散化。

%正弦信号频率分之一是一个周期，采样频率是总的宽度，采样个数是第二个图的幅值。
%当输入信号关于中心对称的时候  IFFT结果为纯的实数,共轭对称，Conjugate symmetry
%When the input signal is symmetrical about the center, the IFFT result is a pure real number.
%实部乘2，虚部抵消
%Real part is multiplied by 2, imaginary part is cancelled


figure();
title("sine");
tiledlayout(4,2);
fs = 100;
X = linspace(0, 2* pi, fs); %生成线性间距向量

Y0 = sin((X + 0) * 1); %f=1, p=0
nexttile;
plot(X,Y0);
title("f=1, p=0");

%幅度和相位和ifft
nexttile;
Yfft = fft(Y0);
Yfft_abs = abs(Yfft);
L_abs = length(Yfft_abs);
%Yfft_abs = abs(Yfft_abs / L_abs); %最大是0-1
Yfft_abs_single = Yfft_abs(1:fix(L_abs/2)+1); %fix == 朝零四舍五入
%Yfft_abs_single(2:end-1) = 2*Yfft_abs_single(2:end-1); %0-1之后除以2了，所以要乘2
FX0=(0:length(Yfft_abs_single)-1);
%m=length(Yfft_abs);
%plot((0:m/2-1)'*fs/m, Yfft_abs(1:m/2));
plot(FX0,Yfft_abs_single);
title("mag");
nexttile;
Yfft_phase = angle(Yfft);
L_phase=length(Yfft_phase);
Yfft_phase_single=Yfft_phase(1:fix(L_phase/2)+1);
%Yfft_phase_single(2:end-1)=2*Yfft_phase_single(2:end-1); %首尾相连
plot(FX0,Yfft_phase_single);
title("phase");
nexttile;
Y_ifft = ifft(Yfft);
plot(X,Y_ifft);
title("ifft");


%改变频率
Y1 = sin((X + 0) * 10); %f=1, p=0
nexttile;
plot(X,Y1);
title("f=10, p=0");

%幅度和相位和ifft
nexttile;
Yfft = fft(Y1);
Yfft_abs = abs(Yfft);
L_abs = length(Yfft_abs);
%Yfft_abs = abs(Yfft_abs / L_abs); %最大是0-1
Yfft_abs_single = Yfft_abs(1:fix(L_abs/2)+1); %fix == 朝零四舍五入
%Yfft_abs_single(2:end-1) = 2*Yfft_abs_single(2:end-1);
FX0=(0:length(Yfft_abs_single)-1);
%m=length(Yfft_abs);
%plot((0:m/2-1)'*fs/m, Yfft_abs(1:m/2));
plot(FX0,Yfft_abs_single);
title("mag");
nexttile;
Yfft_phase = angle(Yfft);
L_phase=length(Yfft_phase);
Yfft_phase_single=Yfft_phase(1:fix(L_phase/2)+1);
%Yfft_phase_single(2:end-1)=2*Yfft_phase_single(2:end-1); %首尾相连
plot(FX0,Yfft_phase_single);
title("phase");
nexttile;
Y_ifft = ifft(Yfft);
plot(X,Y_ifft);
title("ifft");
%% B.3
% a) cos

figure();
title("cos");
tiledlayout(4,2);
fs = 100;
X = linspace(0, 2* pi, fs); %生成线性间距向量

Y0 = cos((X + 0) * 1); %f=1, p=0
nexttile;
plot(X,Y0);
title("f=1, p=0");

%幅度和相位和ifft
nexttile;
Yfft = fft(Y0);
Yfft_abs = abs(Yfft);
L_abs = length(Yfft_abs);
%Yfft_abs = abs(Yfft_abs / L_abs); %最大是0-1
Yfft_abs_single = Yfft_abs(1:fix(L_abs/2)+1); %fix == 朝零四舍五入
%Yfft_abs_single(2:end-1) = 2*Yfft_abs_single(2:end-1); %0-1之后除以2了，所以要乘2
FX0=(0:length(Yfft_abs_single)-1);
%m=length(Yfft_abs);
%plot((0:m/2-1)'*fs/m, Yfft_abs(1:m/2));
plot(FX0,Yfft_abs_single);
title("mag");
nexttile;
Yfft_phase = angle(Yfft);
L_phase=length(Yfft_phase);
Yfft_phase_single=Yfft_phase(1:fix(L_phase/2)+1);
%Yfft_phase_single(2:end-1)=2*Yfft_phase_single(2:end-1); %首尾相连
plot(FX0,Yfft_phase_single);
title("phase");
nexttile;
Y_ifft = ifft(Yfft);
plot(X,Y_ifft);
title("ifft");

%改变频率
Y1 = cos((X + 0) * 10); %f=1, p=0
nexttile;
plot(X,Y1);
title("f=10, p=0");

%幅度和相位和ifft
nexttile;
Yfft = fft(Y1);
Yfft_abs = abs(Yfft);
L_abs = length(Yfft_abs);
%Yfft_abs = abs(Yfft_abs / L_abs); %最大是0-1
Yfft_abs_single = Yfft_abs(1:fix(L_abs/2)+1); %fix == 朝零四舍五入
%Yfft_abs_single(2:end-1) = 2*Yfft_abs_single(2:end-1);
FX0=(0:length(Yfft_abs_single)-1);
%m=length(Yfft_abs);
%plot((0:m/2-1)'*fs/m, Yfft_abs(1:m/2));
plot(FX0,Yfft_abs_single);
title("mag");
nexttile;
Yfft_phase = angle(Yfft);
L_phase=length(Yfft_phase);
Yfft_phase_single=Yfft_phase(1:fix(L_phase/2)+1);
%Yfft_phase_single(2:end-1)=2*Yfft_phase_single(2:end-1); %首尾相连
plot(FX0,Yfft_phase_single);
title("phase");
nexttile;
Y_ifft = ifft(Yfft);
plot(X,Y_ifft);
title("ifft");

% b) Unit Step
figure();
title("unit step");
tiledlayout(2,2);
X = (-1:0.01:1)';
unitstep = X>=0;
nexttile;
plot(X,unitstep);
title("unitstep");
nexttile;
Yfft = fft(unitstep);
Yfft_abs = abs(Yfft);
L_abs = length(Yfft_abs);
%Yfft_abs = abs(Yfft_abs / L_abs); %最大是0-1
Yfft_abs_single = Yfft_abs(1:fix(L_abs/2)+1); %fix == 朝零四舍五入
%Yfft_abs_single(2:end-1) = 2*Yfft_abs_single(2:end-1); %0-1之后除以2了，所以要乘2
FX0=(0:length(Yfft_abs_single)-1);
%m=length(Yfft_abs);
%plot((0:m/2-1)'*fs/m, Yfft_abs(1:m/2));
plot(FX0,Yfft_abs_single);
title("mag");
nexttile;
Yfft_phase = angle(Yfft);
L_phase=length(Yfft_phase);
Yfft_phase_single=Yfft_phase(1:fix(L_phase/2)+1);
%Yfft_phase_single(2:end-1)=2*Yfft_phase_single(2:end-1); %首尾相连
plot(FX0,Yfft_phase_single);
title("phase");
nexttile;
Y_ifft = ifft(Yfft);
plot(X,Y_ifft);
title("ifft");


% c) Impulse
figure();
title("impulse");
tiledlayout(2,2);
%t = (-1:0.001:1)';
impulse = X==0;
nexttile;
plot(X,impulse);
title("impulse");
nexttile;
Yfft = fft(impulse);
Yfft_abs = abs(Yfft);
L_abs = length(Yfft_abs);
%Yfft_abs = abs(Yfft_abs / L_abs); %最大是0-1
Yfft_abs_single = Yfft_abs(1:fix(L_abs/2)+1); %fix == 朝零四舍五入
%Yfft_abs_single(2:end-1) = 2*Yfft_abs_single(2:end-1); %0-1之后除以2了，所以要乘2
FX0=(0:length(Yfft_abs_single)-1);
%m=length(Yfft_abs);
%plot((0:m/2-1)'*fs/m, Yfft_abs(1:m/2));
plot(FX0,Yfft_abs_single);
title("mag");
nexttile;
Yfft_phase = angle(Yfft);
L_phase=length(Yfft_phase);
Yfft_phase_single=Yfft_phase(1:fix(L_phase/2)+1);
%Yfft_phase_single(2:end-1)=2*Yfft_phase_single(2:end-1); %首尾相连
plot(FX0,Yfft_phase_single);
title("phase");
nexttile;
Y_ifft = ifft(Yfft);
plot(X,Y_ifft);
title("ifft");
%% B.4

figure();
tiledlayout(2,2);
phase = pi/2;
f = 1;  %周期等于Π除以f
sample_f = 0.1;
xs = 0:sample_f:2*pi;
ys = 0:sample_f:2*pi;
[X, Y] = meshgrid(xs, ys);
XY = Y*cos(phase) + X*sin(phase);
Z = sin(2*XY*f);
nexttile;
plot3(X, Y, Z);
title("sine wave");
fz=fft2(Z);
fz_abs=abs(fz);
nexttile;
mesh(fz_abs);
title("mag without shift");
fz_shift=fftshift(fz);
fz_shift_abs=abs(fz_shift);
nexttile;
title("mag");
mesh(fz_shift_abs);
nexttile;
ang=angle(fz_shift);
%mesh(ang);   %直接mesh不好看
surf0 = surf(X,Y,ang);
title('phase after shift')
colormap gray
set(surf0, 'edgecolor', 'none') %如果没有这个颜色会变浅或深
view(2) %2D视图
pbaspect([1, 1, 1]);

figure();
tiledlayout(2,2);
phase = pi/2;
f = 2;  %周期等于Π除以f
sample_f = 0.1;
xs = 0:sample_f:2*pi;
ys = 0:sample_f:2*pi;
[X, Y] = meshgrid(xs, ys);
XY = Y*cos(phase) + X*sin(phase);
Z = sin(5*XY);
nexttile;
plot3(X, Y, Z);
title("sine wave");
fz=fft2(Z);
fz_abs=abs(fz);
nexttile;
mesh(fz_abs);
title("mag without shift");
fz_shift=fftshift(fz);
fz_shift_abs=abs(fz_shift);
nexttile;
title("mag");
mesh(fz_shift_abs);
nexttile;
ang=angle(fz_shift);
%mesh(ang);   %直接mesh不好看
surf0 = surf(X,Y,ang);
title('phase after shift')
colormap gray
set(surf0, 'edgecolor', 'none') %如果没有这个颜色会变浅或深
view(2) %2D视图
pbaspect([1, 1, 1]);
%% B.5

% Box Function
figure();
tiledlayout(2,2);
y = 150;
x = 150;
rx = 8;
ry = 8;
[X, Y] = meshgrid(1:1:x, 1:1:y);
Z = zeros([x, y]); %平面是0
x1 = (size(Z, 1) / 2 - rx):(size(Z, 1) / 2 + rx); %size是获取行数和列数
y1 = (size(Z, 2) / 2 - ry):(size(Z, 2) / 2 + ry);
Z(x1, y1) = 1;  %高的地方是1
nexttile;
plot3(X, Y, Z);
title("box function");
fz=fft2(Z);
fz_abs=abs(fz);
nexttile;
mesh(fz_abs);
title("mag without shift");
fz_shift=fftshift(fz);
fz_shift_abs=abs(fz_shift);
nexttile;
mesh(fz_shift_abs);
title("mag");
nexttile;
ang=angle(fz_shift);
mesh(ang);   %直接mesh不好看

%surf0 = surf(X,Y,ang);
%title('phase after shift')
%colormap gray;
%set(surf0, 'edgecolor', 'none') %如果没有这个颜色会变浅或深
%view(2) %2D视图
%pbaspect([1, 1, 1]);


% Gaussian Function
figure();
tiledlayout(2,2);
Z = fspecial('gaussian', 100, 10);
%[X, Y] = meshgrid(1:1:size(Z, 1), 1:1:size(Z, 2));
nexttile;
%plot3(X, Y, Z);
mesh(Z);
title("sine wave");
fz=fft2(Z);
fz_abs=abs(fz);
nexttile;
mesh(fz_abs);
title("mag without shift");
fz_shift=fftshift(fz);
fz_shift_abs=abs(fz_shift);
nexttile;
mesh(fz_shift_abs);
title("mag");
nexttile;
ang=angle(fz_shift);

mesh(ang);   %直接mesh不好看

%surf0 = surf(X,Y,ang);
%title('phase after shift')
%colormap gray;
%set(surf0, 'edgecolor', 'none') %如果没有这个颜色会变浅或深
%view(2) %2D视图
%pbaspect([1, 1, 1]);
%% B.6

figure();
tiledlayout(2,3);
A = imread('C:\Users\Yuhang\Documents\UoS\Image processing\Lab\LenaG.bmp');
nexttile;
imshow(A,[]);
title('original/space domain');
F=fft2(double(A));%将矩阵转化为double型后进行二维傅里叶变换，图像计算中很多处理不能用整型
F_mag = abs(F);
F_phase_without_shift = angle(F);
F_mag_output=log(1+F_mag);%加对数以便于显示图像,为了更好地显示细节，不进行log变换的话灰度的动态范围被压缩

nexttile;
imshow(F_mag_output,[]);%根据 C 中像素值的范围缩放显示.使用 [min（C（:)）max（C（:)）] 作为显示范围。 imshow 将 C 中的最小值显示为黑色，将最大值显示为白色。使用imshow(A,[])，即可把图像矩阵A显示为正常的灰度图像。本来A是0-1，把double拉伸到[0 255]
title('spectrum magnitude without shift');
nexttile;
imshow(F_phase_without_shift,[]);
title('phase without shift');

F_shift=fftshift(F);%fftshift将傅里叶变换的零频率成分移到频谱中心，因为fft2变换中，信号的零频率成分在信号左上角。
F_mag_shift=log(1+abs(F_shift));
F_phase_shift=angle(F_shift);
nexttile;
imshow(F_mag_shift,[]);
title('spectrum magnitude with shift');
nexttile;
imshow(F_phase_shift,[]);
title('phase with shift');
%% B.7

figure();
tiledlayout(2,3);
A = imread('C:\Users\Yuhang\Documents\UoS\Image processing\Lab\LenaG.bmp');
nexttile;
imshow(A,[]);
title('original/space domain');
F=fft2(double(A));%将矩阵转化为double型后进行二维傅里叶变换，图像计算中很多处理不能用整型
F_mag = abs(F);
F_phase_without_shift = angle(F);
F_mag_output=log(1+F_mag);%加对数以便于显示图像,为了更好地显示细节，不进行log变换的话灰度的动态范围被压缩

nexttile;
imshow(F_mag_output,[]);%[]和mat2gray一个意思，根据 C 中像素值的范围缩放显示.使用 [min（C（:)）max（C（:)）] 作为显示范围。 imshow 将 C 中的最小值显示为黑色，将最大值显示为白色。使用imshow(A,[])，即可把图像矩阵A显示为正常的灰度图像。本来A是0-1，把double拉伸到[0 255]
title('spectrum magnitude without shift');
nexttile;
imshow(F_phase_without_shift,[]);
title('phase without shift');

F_shift=fftshift(F);%fftshift将傅里叶变换的零频率成分移到频谱中心，因为fft2变换中，信号的零频率成分在信号左上角。
F_mag_shift=log(1+abs(F_shift));
F_phase_shift=angle(F_shift);
nexttile;
imshow(F_mag_shift,[]);
title('spectrum magnitude with shift');
nexttile;
imshow(F_phase_shift,[]);
title('phase with shift');

nexttile;
F_ifft = ifft2(F);
imshow(F_ifft,[]);
title("image after ifft");

%虚部代表相位延后90度即时间延后四分之一个周期

%实部乘2，虚部抵消
%Real part is multiplied by 2, imaginary part is cancelled
%% B.8

figure();
tiledlayout(2,2);
A = imread('C:\Users\Yuhang\Documents\UoS\Image processing\Lab\LenaG.bmp');
nexttile;
imshow(A,[]);
title('original/space domain');
F=fft2(double(A));%将矩阵转化为double型后进行二维傅里叶变换，图像计算中很多处理不能用整型
F_mag = abs(F);
F_phase = angle(F);

nexttile;
F_ifft = ifft2(F);
imshow(F_ifft,[]);
title("image after ifft");

%New part
phase0= F_phase .* 0;
F_phase0=F_mag .* cos(phase0) + F_mag .* sin(phase0) .* 1i;
F_phase0_ifft2=log((abs(ifft2(F_phase0))));
nexttile;
imshow(F_phase0_ifft2,[]);
title("all phase == 0, without shift");

nexttile;
%imshow(mat2gray(log(fftshift(bmp_img_inv_zero_phase))))
imshow(log(fftshift(ifft2(F_mag))),[]); % []和mat2gray:归一化，让每个元素的值都在0和1之间
title("2nd way but with shift");


% and the phase determines the position information.


%% B.9

figure();
tiledlayout(2,2);
A = imread('C:\Users\Yuhang\Documents\UoS\Image processing\Lab\LenaG.bmp');
nexttile;
imshow(A,[]);
title('original/space domain');
F=fft2(double(A));%将矩阵转化为double型后进行二维傅里叶变换，图像计算中很多处理不能用整型
F_mag = abs(F);
F_phase = angle(F);

nexttile;
F_ifft = ifft2(F);
imshow(F_ifft,[]);
title("image after ifft");

%New part
mag1= F_mag ./ F_mag;
%mag1 = max(abs(F), [], 'all'); %max查找最大值
F_mag1=mag1 .* cos(F_phase) + mag1 .* sin(F_phase) .* 1i;
F_mag1_ifft2=((ifft2(F_mag1)));
nexttile;
imshow(mat2gray(real(F_mag1_ifft2))); %mat2gray必须对实数操作，而且，必须显示实数，有虚数可能是因为误差。
title("unit mag");

%The amplitude in the space determines the intensity information,so the contour can be seen.
%% B.10;

figure();
tiledlayout(2,3);
A = imread('C:\Users\Yuhang\Documents\UoS\Image processing\Lab\LenaG.bmp');

%A=rot90(A); %应该用哪一个呢？
A=A';

nexttile;
imshow(A,[]);
title('original/space domain');
F=fft2(double(A));%将矩阵转化为double型后进行二维傅里叶变换，图像计算中很多处理不能用整型
F_mag = abs(F);
F_phase_without_shift = angle(F);
F_mag_output=log(1+F_mag);%加对数以便于显示图像,为了更好地显示细节，不进行log变换的话灰度的动态范围被压缩

nexttile;
imshow(F_mag_output,[]);%根据 C 中像素值的范围缩放显示.使用 [min（C（:)）max（C（:)）] 作为显示范围。 imshow 将 C 中的最小值显示为黑色，将最大值显示为白色。使用imshow(A,[])，即可把图像矩阵A显示为正常的灰度图像。本来A是0-1，把double拉伸到[0 255]
title('spectrum magnitude without shift');
nexttile;
imshow(F_phase_without_shift,[]);
title('phase without shift');

F_shift=fftshift(F);%fftshift将傅里叶变换的零频率成分移到频谱中心，因为fft2变换中，信号的零频率成分在信号左上角。
F_mag_shift=log(1+abs(F_shift));
F_phase_shift=angle(F_shift);
nexttile;
imshow(F_mag_shift,[]);
title('spectrum magnitude with shift');
nexttile;
imshow(F_phase_shift,[]);
title('phase with shift');

%When the original image is rotated, the Fourier transform will also rotate. 
%The decomposition to the spatial frequency reflects the orientation of the features in the image, so the Fourier transform has orientation dependence