%% HW4 Program Assignment
% By: 105060012 張育菘
%
%% P1
% Let x\left\lbrack n\right\rbrack ={n\left(0\ldotp 9\right)}^n u\left\lbrack n\right\rbrack,
%
% (a) Determine the DTFT \widetilde{X} \left(e^{j\omega } \right) of x\left\lbrack n\right\rbrack. Please write your calculations and answer on your .mlx file.
% Ans.
% a^n u\left\lbrack n\right\rbrack \leftrightarrow \frac{1}{1-{\mathrm{az}}^{-1} },\mathrm{nx}\left\lbrack n\right\rbrack \leftrightarrow -z\;\frac{\mathrm{dX}\left(z\right)}{\mathrm{dz}},z=e^{j\omega } \Longrightarrow {\mathrm{na}}^n u\left\lbrack n\right\rbrack \leftrightarrow \frac{{\mathrm{az}}^{-1\;} }{{\left(1-{\mathrm{az}}^{-1} \right)}^2 }=\frac{ae^{-j\omega } }{{\left(1-ae^{-j\omega } \right)}^2 }
% \Longrightarrow {n\left(0\ldotp 9\right)}^n u\left\lbrack n\right\rbrack \leftrightarrow \frac{0\ldotp 9e^{-j\omega } }{{\left(1-0\ldotp 9e^{-j\omega } \right)}^2 }
%
% (b) Choose first N = 20 samples of x[n] and compute the approximate DTFT \tilde{X_N } \left(e^{j\omega } \right) using the fft function. Plot magnitudes of \tilde{X} \left(e^{j\omega } \right) and \tilde{X_N } \left(e^{j\omega } \right) in one plot and compare your results.

close all; clear;
fprintf('1(b)\n');

n = linspace(0, 1023, 1024);

x = n.*((0.9).^n);
X_ft = fft(x); X_ft = fftshift(X_ft);
om = linspace(-pi, pi, 1024); 
X10_ft = fft(x, 10); X10_ft = fftshift(X10_ft);
om_b = linspace(-pi, pi, 10);

figure; 
plot(om_b/pi,abs(X10_ft)); hold on;
plot(om/pi, abs(X_ft)); hold off; grid on;

title('Magnitude Response of DTFT with N = 10');
legend('X_1_0(e^j^\omega)','X(e^j^\omega)')
ylabel('Magnitude');
xlabel('Time(pi)');

%
% (c) Repeat part (b) using N = 50.

fprintf('1(c)\n');

X50_ft = fft(x, 50); X50_ft = fftshift(X50_ft);
om_c = linspace(-pi, pi, 50);

figure; 
plot(om_c/pi,abs(X50_ft)); hold on;
plot(om/pi, abs(X_ft)); hold off; grid on;

title('Magnitude Response of DTFT with N = 50');
legend('X_5_0(e^j^\omega)','X(e^j^\omega)')
ylabel('Magnitude');
xlabel('Time(pi)');

%
% (d) Repeat part (b) using N = 100.

fprintf('1(d)\n');

X100_ft = fft(x, 100); X100_ft = fftshift(X100_ft);
om_d = linspace(-pi, pi, 100);

figure; 
plot(om_d/pi,abs(X100_ft)); hold on;
plot(om/pi, abs(X_ft)); hold off; grid on;

title('Magnitude Response of DTFT with N = 100');
legend('X_1_0_0(e^j^\omega)','X(e^j^\omega)')
ylabel('Magnitude');
xlabel('Time(pi)');

%
%% P2
% Let x\left\lbrack n\right\rbrack =x_1 \left\lbrack n\right\rbrack +{\mathrm{jx}}_2 \left\lbrack n\right\rbrack where sequences x_1 \left\lbrack n\right\rbrack and x_2 \left\lbrack n\right\rbrack are real-valued.
%
% (a)Show that X_1 \left\lbrack k\right\rbrack =X^{\mathrm{cce}} \left\lbrack k\right\rbrack and {\mathrm{jX}}_2 \left\lbrack k\right\rbrack =X^{\mathrm{cco}} \left\lbrack k\right\rbrack. Please write your calculations and answer on your .mlx file.
% Ans.
% X^{\mathrm{cce}} \left\lbrack k\right\rbrack =\frac{1}{2}\left(X\left\lbrack k\right\rbrack +X^* \left\lbrack {\left\langle -k\right\rangle }_N \right\rbrack \right),X^{\mathrm{cco}} \left\lbrack k\right\rbrack =\frac{1}{2}\left(X\left\lbrack k\right\rbrack -X^* \left\lbrack {\left\langle -k\right\rangle }_N \right\rbrack \right),\;x^* \left\lbrack n\right\rbrack \leftrightarrow X^* \left\lbrack {\left\langle -k\right\rangle }_N \right\rbrack
% x\left\lbrack n\right\rbrack =x_1^{\textrm{ce}} \left\lbrack n\right\rbrack +x_1^{\textrm{co}} \left\lbrack n\right\rbrack +{\textrm{jx}}_2^{\textrm{ce}} \left\lbrack n\right\rbrack +{\textrm{jx}}_2^{\textrm{co}} \left\lbrack n\right\rbrack \leftrightarrow X\left\lbrack k\right\rbrack =X_1^{\textrm{ce}} \left\lbrack k\right\rbrack +{\textrm{jX}}_1^{\textrm{co}} \left\lbrack k\right\rbrack +{\textrm{jX}}_2^{\textrm{ce}} \left\lbrack k\right\rbrack +X_2^{\textrm{co}} \left\lbrack k\right\rbrack
% x^* \left\lbrack n\right\rbrack =x_1^{\textrm{ce}} \left\lbrack n\right\rbrack +x_1^{\textrm{co}} \left\lbrack n\right\rbrack -{\textrm{jx}}_2^{\textrm{ce}} \left\lbrack n\right\rbrack -{\textrm{jx}}_2^{\textrm{co}} \left\lbrack n\right\rbrack \leftrightarrow X^* \left\lbrack {\left\langle -k\right\rangle }_N \right\rbrack =X_1^{\textrm{ce}} \left\lbrack k\right\rbrack +{\textrm{jX}}_1^{\textrm{co}} \left\lbrack k\right\rbrack -{\textrm{jX}}_2^{\textrm{ce}} \left\lbrack k\right\rbrack -X_2^{\textrm{co}} \left\lbrack k\right\rbrack
% \Longrightarrow X_1 \left\lbrack k\right\rbrack =X_1^{\mathrm{ce}} \left\lbrack k\right\rbrack +{\mathrm{jX}}_1^{\mathrm{co}} \left\lbrack k\right\rbrack =\frac{1}{2}\left(X\left\lbrack k\right\rbrack +X^* \left\lbrack {\left\langle -k\right\rangle }_N \right\rbrack \right)=X^{\mathrm{cce}} \left\lbrack k\right\rbrack
% \Longrightarrow {\mathrm{jX}}_2 \left\lbrack k\right\rbrack ={\textrm{jX}}_2^{\textrm{ce}} \left\lbrack k\right\rbrack +X_2^{\textrm{co}} \left\lbrack k\right\rbrack =\frac{1}{2}\left(X\left\lbrack k\right\rbrack -X^* \left\lbrack {\left\langle -k\right\rangle }_N \right\rbrack \right)=X^{\textrm{cco}} \left\lbrack k\right\rbrack
%
% (b) Write a MATLAB function [X1,X2] = tworealDFTs(x1,x2) that implements the results in part (a).

close all; clear;
fprintf('2(b)\n');

open tworealDFTs.m;

%
% (c) Verify your function on the following two sequences: x_1 \left\lbrack n\right\rbrack =0\ldotp 9^n ,x_2 \left\lbrack n\right\rbrack =\left(1-0\ldotp 8^n \right);0\le n\le 49

fprintf('2(c)\n');

N = 50;
n = linspace(0, N-1, N);
x1 = ((0.9).^n);
x2 = (1-((0.8).^n));

[X1_ft,X2_ft] = tworealDFTs(x1,x2)
X1_ft = fftshift(X1_ft); X2_ft = fftshift(X2_ft);
om = linspace(-pi, pi, N);

X1_fft = fft(x1); X1_fft = fftshift(X1_fft);
X2_fft = fft(x2); X2_fft = fftshift(X2_fft);

figure; 
plot(om/pi,abs(X1_fft)); hold on;
plot(om/pi,abs(X1_ft)); hold off; grid on;
title('Magnitude Response of DTFT of X1(e^j^\omega)');
legend('by fft','by tworealDFTs');
ylabel('Magnitude');
xlabel('Time(pi)');

figure; 
plot(om/pi,abs(X2_fft)); hold on;
plot(om/pi,abs(X2_ft)); hold off; grid on;
title('Magnitude Response of DTFT X2(e^j^\omega)');
legend('by fft','by tworealDFTs');
ylabel('Magnitude');
xlabel('Time(pi)');

%
%% P3
% Let x_1 \left\lbrack n\right\rbrack =\left\lbrace 1_{n=0} ,2,3,4,5\right\rbracebe a 5-point sequence and let x_2 \left\lbrack n\right\rbrack =\left\lbrace 2_{n=0} ,-1,1,-1\right\rbrace be a 4-point sequence.
%
% (a) Determine x_1 \left\lbrack n\right\rbrack \odot x_2 \left\lbrack n\right\rbrack using hand calculations. Please write your calculations and answer on your .mlx file.
% Ans.
% If length(x_1 \left\lbrack n\right\rbrack) \not= length(x_2 \left\lbrack n\right\rbrack) \Longrightarrowzero padding
% x_1 \left\lbrack n\right\rbrack =\left\lbrace 1_{n=0} ,2,3,4,5\right\rbrace ,x_2 \left\lbrack n\right\rbrack =\left\lbrace 2_{n=0} ,-1,1,-1,\;0\right\rbrace
% y\left\lbrack n\right\rbrack =x_1 \left\lbrack n\right\rbrack \odot x_2 \left\lbrack n\right\rbrack =\sum_{m\;=\;0}^4 x_1 \left\lbrack m\right\rbrack x_2 \left\lbrack {\left\langle n-m\right\rangle }_5 \right\rbrack
% \begin{array}{l}
% \Longrightarrow y\left\lbrack 0\right\rbrack =x_1 \left\lbrack 0\right\rbrack x_2 \left\lbrack 0\right\rbrack +x_1 \left\lbrack 1\right\rbrack x_2 \left\lbrack -1\right\rbrack +x_1 \left\lbrack 2\right\rbrack x_2 \left\lbrack -2\right\rbrack +x_1 \left\lbrack 3\right\rbrack x_2 \left\lbrack -3\right\rbrack +x_1 \left\lbrack 4\right\rbrack x_2 \left\lbrack -4\right\rbrack =2+0-3+4-5=-2\\
% \;\;\;\;\;\;\;\;\;y\left\lbrack 1\right\rbrack =x_1 \left\lbrack 0\right\rbrack x_2 \left\lbrack 1\right\rbrack +x_1 \left\lbrack 1\right\rbrack x_2 \left\lbrack 0\right\rbrack +x_1 \left\lbrack 2\right\rbrack x_2 \left\lbrack -1\right\rbrack +x_1 \left\lbrack 3\right\rbrack x_2 \left\lbrack -2\right\rbrack +x_1 \left\lbrack 4\right\rbrack x_2 \left\lbrack -3\right\rbrack =-1+4+0-4+5=4\\
% \;\;\;\;\;\;\;\;\;y\left\lbrack 2\right\rbrack =x_1 \left\lbrack 0\right\rbrack x_2 \left\lbrack 2\right\rbrack +x_1 \left\lbrack 1\right\rbrack x_2 \left\lbrack 1\right\rbrack +x_1 \left\lbrack 2\right\rbrack x_2 \left\lbrack 0\right\rbrack +x_1 \left\lbrack 3\right\rbrack x_2 \left\lbrack -1\right\rbrack +x_1 \left\lbrack 4\right\rbrack x_2 \left\lbrack -2\right\rbrack =1-2+6+0-5=0\\
% \;\;\;\;\;\;\;\;\;y\left\lbrack 3\right\rbrack =x_1 \left\lbrack 0\right\rbrack x_2 \left\lbrack 3\right\rbrack +x_1 \left\lbrack 1\right\rbrack x_2 \left\lbrack 2\right\rbrack +x_1 \left\lbrack 2\right\rbrack x_2 \left\lbrack 1\right\rbrack +x_1 \left\lbrack 3\right\rbrack x_2 \left\lbrack 0\right\rbrack +x_1 \left\lbrack 4\right\rbrack x_2 \left\lbrack -1\right\rbrack =-1+2-3+8+0=6\\
% \;\;\;\;\;\;\;\;\;y\left\lbrack 4\right\rbrack =x_1 \left\lbrack 0\right\rbrack x_2 \left\lbrack 4\right\rbrack +x_1 \left\lbrack 1\right\rbrack x_2 \left\lbrack 3\right\rbrack +x_1 \left\lbrack 2\right\rbrack x_2 \left\lbrack 2\right\rbrack +x_1 \left\lbrack 3\right\rbrack x_2 \left\lbrack 1\right\rbrack +x_1 \left\lbrack 4\right\rbrack x_2 \left\lbrack 0\right\rbrack =0-2+3-4+10=7\\
% \Longrightarrow y\left\lbrack n\right\rbrack =\left\lbrace {-2}_{n=0} ,4,0,6,7\right\rbrace
% \end{array}
%
% (b) Verify your calculations in (a) using the circonv function.

close all; clear;
fprintf('3(b)\n');

x1 = [1 2 3 4 5]'; 
x2 = [2 -1 1 -1]'; 
x2n = [2 -1 1 -1 0]'; % zero padding
y = circonv(x1, x2n);
y = y'

% We can find that the results of 3(a) and 3(b) are the same.
%
% (c) Verify your calculations in (a) by computing the DFTs and IDFT.

fprintf('3(c)\n');
ydft = ifft(fft(x1).*fft(x2n));
ydft = ydft'

% We can find that the results of 3(a) and 3(c) are the same.
%
%% P4
% Let x_1 \left\lbrack n\right\rbrack be an N_1-point and x_2 \left\lbrack n\right\rbrack be an N_2-point sequence. Let N ≥ max(N1,N2). Their N-point circular convolution is shown to be equal to the aliased version of their linear convolution in (7.209) in Program Assignment 3. This result can be used to compute the circular convolution via the linear convolution.
% x_3 \left\lbrack n\right\rbrack =x_1 \left\lbrack n\right\rbrack *x_2 \left\lbrack n\right\rbrack ,x_4 \left\lbrack n\right\rbrack =x_1 \left\lbrack n\right\rbrack \odot x_2 \left\lbrack n\right\rbrack
% x_4 \left\lbrack n\right\rbrack =\sum_{l=-\infty }^{\infty } x_3 \left\lbrack n+\mathrm{lN}\right\rbrack \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left(7\ldotp 209\right)
%
% (a) Develop a MATLAB function y = lin2circonv(x,h) that implements this approach.

close all; clear;
fprintf('4(a)\n');
open lin2circonv.m;

%
% (b) For x\left\lbrack n\right\rbrack =\left\lbrace 1_{n=0} ,2,3,4\right\rbrace and h\left\lbrack n\right\rbrack =\left\lbrace 1_{n=0} ,-1,1,-1\right\rbrace determine their 4-point circular convolution using the lin2circonv function and verify using the circonv function.

fprintf('4(b)\n');

x = [1 2 3 4]'; h = [1 -1 1 -1]';
hx_circonv = circonv(h, x);
hx_circonv = hx_circonv'
hx_lin2circonv = lin2circonv(h', x')

% We can find that the results by the function circonv and lin2circonv are the same.
%
%% P5
% Let a 2D filter impulse response h[m, n] be given by
% h\left\lbrack m,n\right\rbrack =\left\lbrace \begin{array}{cc}
% \frac{1}{2\pi \sigma^2 }e^{-\frac{m^2 +n^2 }{2\sigma^2 }}  & ,-128\le m,n\le 127\\
% 0 & ,\textrm{otherwise}
% \end{array}\right.
% where σ is a parameter. For this problem use the “Lena” image.
%
% (a) For σ = 4, determine h[m, n] and compute its 2D-DFT H[k, l] via the fft2 function taking care of shifting the origin of the array from the middle to the beginning (using the ifftshift function). Show the log-magnitude of H[k, l] as an image.

close all; clear;
fprintf('5(a)\n');

sigma = 4;
h = zeros(256, 256);

for m = -128:1:127
    for n = -128:1:127
        h(m+129, n+129) = ((2*pi*(sigma)^2)^-1).*exp(-1*((m^2)+(n^2))./(2*(sigma)^2));
    end
end

H = ifftshift(h);
H = fft2(H)
H_abs = abs(H);
H_log = log(H_abs);

figure; imshow(H_log,[]); title('Filter with \sigma = 4');
figure; imagesc(H_log);

%
% (b) Process the “Lena” image in the frequency domain using the above H[k, l]. This will involve taking 2D-DFT of the image, multiplying the two DFTs and then taking the inverse of the product. Comment on the visual quality of the resulting filtered image.

fprintf('5(b)\n');

sigma = 4;
h = zeros(256, 256);

for m = -128:1:127
    for n = -128:1:127
        h(m+129, n+129) = ((2*pi*(sigma)^2)^-1).*exp(-1*((m^2)+(n^2))./(2*(sigma)^2));
    end
end

H = ifftshift(h); H = fft2(H);

lena = imread('lena.jpg');
lena_fft2 = fft2(lena);
lena_H4 = lena_fft2.*H;
lena_h4_ifft = ifft2(lena_H4);
figure; imshow(abs(lena)); title('Original version');
figure; imshow(abs(lena_h4_ifft),[]); title('Afer filtering (\sigma = 4)');

% We find that the image after filtering is more blurred than the original image.
%
% (c) Repeat (a) and (b) for σ = 32 and comment on the resulting filtered image as well as the difference between the two filtered images.

fprintf('5(c)\n');

sigma = 32;
h = zeros(256, 256);

for m = -128:1:127
    for n = -128:1:127
        h(m+129, n+129) = ((2*pi*(sigma)^2)^-1).*exp(-1*((m^2)+(n^2))./(2*(sigma)^2));
    end
end

H = ifftshift(h); H = fft2(H);
H_abs = abs(H); H_log = log(H_abs);
figure; imshow(H_log,[]); title('Filter with \sigma = 32');

lena = imread('lena.jpg');
lena_fft2 = fft2(lena);
lena_H32 = lena_fft2.*H;
lena_h32_ifft = ifft2(lena_H32);
% figure; imshow(abs(lena)); title('Original version');
figure; imshow(abs(lena_h32_ifft),[]); title('Afer filtering (\sigma = 32)');

% 比較 sigma = 4 與 sigma = 32 的 filter，我們可以發現 filter 從四角是白色變成像是兩片為透明的黑色片沒有對準疊在一起的圖形。
% Compared with the previous image with sigma 4, and this image with sigma 32, we can find that the image with sigma 32 is more blurred than the previous image.
%
% (d) The filtered image in part (c) also suffers from an additional distortion due to a spatial-domain aliasing effect in the circular convolution. To eliminate this artifact, consider both the image and the filter h[m, n] as 512 × 512 size images using zero-padding in each dimension. Now perform the frequency-domain filtering and comment on the resulting filtered image.

fprintf('5(d)\n');

sigma = 32;
h = zeros(256, 256);
for m = -128:1:127
    for n = -128:1:127
        h(m+129, n+129) = ((2*pi*(sigma)^2)^-1).*exp(-1*((m^2)+(n^2))./(2*(sigma)^2));
    end
end

h_5d = padarray(h, [128,128]);

H_5d = ifftshift(h_5d);
H_5d = fft2(H_5d)
H_5d_abs = abs(H_5d); H_5d_log = log(H_5d_abs);
figure; imshow(H_5d_log,[]); title('Filter with zero-padding and \sigma = 32');

lena = imread('lena.jpg');
lena_d = padarray(lena, [128 128]);
lena_fft2 = fft2(lena_d);
lena_H32d = lena_fft2.*H_5d;
lena_h32d_ifft = ifft2(lena_H32d);
% figure; imshow(abs(lena_d)); title('Original version');
figure; imshow(abs(lena_h32d_ifft),[]); title('Afer zero-padding filter (\sigma = 32)');

% 我們發現經過 zero-padding 處理過的圖形較沒有經過 zero-padding 處理的圖形來的模糊，甚而言之，有經過 zero-padding 處理過的圖形有邊界扭曲的情形發生。
%
% (e) Repeat part (b) for σ = 4 but now using the frequency response 1 − H[k, l] for the filtering. Compare the resulting filtered image with that in (b).

fprintf('5(e)\n');

sigma = 4;
h = zeros(256, 256);
for m = -128:1:127
    for n = -128:1:127
        h(m+129, n+129) = ((2*pi*(sigma)^2)^-1).*exp(-1*((m^2)+(n^2))./(2*(sigma)^2));
    end
end

H = ifftshift(h);
H = fft2(H);
H_5e = 1 - H;
H_5e_abs = abs(H_5e); H_5e_log = log(H_5e_abs);
figure; imshow(H_5e_log,[]); title('1 - H Filter with \sigma = 4');

lena = imread('lena.jpg');
lena_fft2 = fft2(lena);
lena_H4e = lena_fft2.*H_5e;
lena_h4e_ifft = ifft2(lena_H4e);
% figure; imshow(abs(lena)); title('Original version');
figure; imshow(abs(lena_h4e_ifft),[]); title('Afer (1 - H) filter (\sigma = 4)');

% 從 5_e 的結果可發現 5_d 的圖形的人物輪廓比較模糊，但是同時圖片比較白。