% Table of Contents
% Filtering in the Frequency Domain
%             Low Pass Filter
%     Frequency Space Filtering Methods
%         Low Pass Filtering with the Fourier Transform:
%             Ideal 2D Low-Pass Filter
%             The Not So Ideal Low-Pass Filter?
%         MATLAB Ideal Low Pass Filter Example 1
%     MATLAB Ideal Low-Pass Filter Example 2
%         Low-Pass Butterworth Filter
%             Visualising the 1D Low-Pass Butterworth Filter
%         MATLAB Butterworth Low Pass Filter Example 1
%         MATLAB Low-Pass Butterworth Filter Example 2
%         Low Pass Filtering Noisy Images
%     Other Filters
%         High Pass Filtering
% Convolution
%     1D Convolution
%         1D Convolution Example in the Real Domain
%     The Fourier Transform and Convolution
%             Computing Convolutions with the Fourier Transform
%     Image Deblurring Deconvolution Example
%     A Practical Audio Example of Convolution: Reverb
%     Reverb Simulation
%     Schroeder's Reverberator
%         Schroeder's reverberator example (one of a few variations):
%         A simple Schroeder reverberator,  schroeder1.m:
%     Convolution Reverb Implementation
%     Room Impulse Responses
%     MATLAB Convolution Reverb
% Further Reading
%     Supporting functions
%%
%% Filtering in the Frequency Domain
%% Low Pass Filter
%
% Example uses of a Low Pass Filter: To remove Audio Hiss, 'Salt and Pepper' noise in images,
%
% Noise filtering: The idea with noise Filtering is to reduce various spurious effects of a local nature in the image, caused perhaps by
% noise in the acquisition system,
% arising as a result of transmission of the data, for example from a space probe utilising a low-power transmitter.
%
%%
%
%% Frequency Space Filtering Methods
%
% Low Pass Filtering — Aim to Remove Noise
% Noise = High Frequencies:
% In audio data many spurious peaks in over a short timescale.
% In an image means there are many rapid transitions (over a short distance) in intensity from high to low and back again or vice versa, as faulty pixels are encountered.
% Not all high frequency data noise though!
% Therefore,  noise will contribute heavily to the high frequency components of the signal when it is analysed in Fourier space.
% Thus if we reduce the high frequency components — Low-Pass Filter should (if tuned properly) reduce the amount of noise in the data.
%
%% Low Pass Filtering with the Fourier Transform:
%  We filter in Fourier space by computing G(u,v)=H(u,v) F(u,v) where:
% F(u,v) is the Fourier transform of the original image,
% H(u,v) is a filter function, designed to reduce high frequencies, and
% G(u,v) is the Fourier transform of the improved image.
% Inverse Fourier transform G(u,v) to get g(x,y) our improved image
% We need to design or compute the filter function, H(u,v). Two ways:
% If we know h(x,y) or have a discrete sample of h(x,y) can compute its Fourier Transform
% Can simply design/craft simple filters in Frequency Space
% The simplest sort of filter to use is an ideal low-pass filter, which in one dimension appears as:
%
%%
%
% This is a top-hat function (recall realtion to Sinc function above), H(u), which is 1 for u between 0 and u_0, the cut-off frequency, and zero elsewhere.
% So all frequency space information above u_0 is discarded, and all information below u_0 is kept.
% A very simple computational process.
%
%% Ideal 2D Low-Pass Filter
%
% The two dimensional version of this is the Low-Pass Filter:
% H(u,v) = 1 & \mbox{if $\sqrt{u^2+v^2} \leq w_0$}\\ 0 & \mbox{otherwise,} \end{array} \right. $  if \sqrt{u^2+v^2} \leq w_0
%        0   otherwise
% where w_0 is now the cut-off frequency for both dimensions.
% Thus, all frequencies inside a radius w_0 are kept, and all others discarded.
%
%
%% The Not So Ideal Low-Pass Filter?
%
% In practice, the ideal Low-Pass Filter is not so ideal.
% The problem with this filter is that as well as noise there may be useful high frequency content:
% In audio: plenty of other high frequency content: high pitches, rustles, scrapes, wind, mechanical noises, cymbal crashes etc.
% In images: edges (places of rapid transition from light to dark) also significantly contribute to the high frequency components.
% Choosing the most appropriate cut-off frequency is not so easy
% Similar problem to choosing a threshold in image thresholding.
% So what happens if you set the wrong value for the cut-off frequency?
% If you choose the wrong cut-off frequency an ideal low-pass filter will tend to  blur the data:
% High audio frequencies become muffled
% Edges in images become blurred.
% The lower the cut-off frequency is made, the more pronounced this effect becomes in useful data content.
%
%% MATLAB Ideal Low Pass Filter Example 1
% lowpass.m:

% Create a white box on a black background image
M= 256; N = 256;
image = zeros(M,N);
box = ones(64,64);
%box at centre
image(97:160,97:160) = box;

% Show Image

figure
imshow(image)
title('Original Image');

% compute fft and display its spectra

F=fft2(double(image));
figure
imagesc((abs(fftshift(F))/(M*N)));
colormap(jet);
axis off;
title('FFT Image Spectra');

%compute Ideal Low Pass Filter

u0 = 20; % set cut off frequency


u=0:(M-1);
v=0:(N-1);
idx=find(u>M/2);
u(idx)=u(idx)-M;
idy=find(v>N/2);
v(idy)=v(idy)-N;
[V,U]=meshgrid(v,u);
D=sqrt(U.^2+V.^2);
H=double(D<=u0);

% display
figure
imshow(fftshift(H))
title('Ideal Low Pass Filter');

% Apply filter and do inverse FFT
G=H.*F;
g=real(ifft2(double(G)));

% Show Result

figure
imshow(g);
title('Low Pass Filtered Image');

%%% Demo now to show noisy filtering properties of the low pass filter

% now addnoise to the image

image_noisy = imnoise(image,'gaussian');

% show noisy im
figure
imshow(image_noisy)
title('Image with Added Noise');


%%%% compute FFT of noisy image

F = fft2(double(image_noisy));



% Filter using low pass filter H as above

G_lowpass=H.*F;

% convert back to image spce with iFFT

g_lowpass = real(ifft2(double(G_lowpass)));

figure
imshow(g_lowpass);
title('Low Pass Filtered Noisy Image');

%
%% MATLAB Ideal Low-Pass Filter Example 2
% Same example as above but a different image with a differently tuned cut off frequencies: compare the difference in output.
% lowpass2.m:

% read in MATLAB demo text image
image = imread('text.png');
[M, N] = size(image);


% Show Image

figure;
imshow(image);
title('Original Image');

% compute fft and display its spectra

F=fft2(double(image));
figure;
imagesc((abs(fftshift(F))/(M*N)));
colormap(jet);
axis off;
title('FFT Image Spectra');

%compute Ideal Low Pass Filter

u0 = 50; % set cut off frequency


u=0:(M-1);
v=0:(N-1);
idx=find(u>M/2);
u(idx)=u(idx)-M;
idy=find(v>N/2);
v(idy)=v(idy)-N;
[V,U]=meshgrid(v,u);
D=sqrt(U.^2+V.^2);
H=double(D<=u0);

% display
figure;
imshow(fftshift(H));
title('Ideal Low Pass Filter');

% Apply filter and do inverse FFT
G=H.*F;
g=real(ifft2(double(G)));

% Show Result

figure;
imshow(g);
title('Low Pass Filtered Image');
%%% Demo now to show noisy filtering properties of the low pass filter

% now addnoise to the image

image_noisy = imnoise(double(image),'gaussian');

% show noisy im
figure
imshow(image_noisy)
title('Image with Added Noise');


%%%% compute FFT of noisy image

F = fft2(double(image_noisy));



% Filter using low pass filter H as above

G_lowpass=H.*F;

% convert back to image spce with iFFT

g_lowpass = real(ifft2(double(G_lowpass)));

figure
imshow(g_lowpass);
title('Low Pass Filtered Noisy Image');



%%%%% Now set Higher cut off frequency and repeat all above

%compute Ideal Low Pass Filter

u0 = 80; % set cut off frequency


u=0:(M-1);
v=0:(N-1);
idx=find(u>M/2);
u(idx)=u(idx)-M;
idy=find(v>N/2);
v(idy)=v(idy)-N;
[V,U]=meshgrid(v,u);
D=sqrt(U.^2+V.^2);
H=double(D<=u0);

% display
figure;
imshow(fftshift(H));
title('High Cut-off Frequency Low Pass Filter');

% Apply filter and do inverse FFT


% Filter using low pass filter H as above

G_lowpass=H.*F;

% convert back to image spce with iFFT 

g_lowpass = real(ifft2(double(G_lowpass)));

figure
imshow(g_lowpass);
title('High Cut-off Frequency Low Pass Filtered Image');

%
%% Low-Pass Butterworth Filter
% Another popular (and general) filter is the Butterworth low pass filter.
% We introduced the Butterworth Filter with IIR/FIR Filters (Temporal Domain Filtering).
% Let's now explain this it in more detail.
% Much easier to visualise in Frequency space
% In the 2D case, H(u,v) takes the form
%  H(u,v)=\frac{1}{1+\left[(u^2+v^2)/w_0^2\right]^n},
% where n is called the order of the filter.
%
%
%% Visualising the 1D Low-Pass Butterworth Filter
% This keeps some of the high frequency information, as illustrated by the second order one dimensional Butterworth filter:
%
% Consequently reduces the blurring.
% Blurring the filter — Butterworth is essentially a smoothed top hat functions — reduces blurring by the filter.
% Visualising the 2D Low-Pass Butterworth Filter The 2D second order Butterworth filter looks like this:
%
%
% Note this is blurred circle — blurring of the ideal 2D Low-Pass Filter.
%
%% MATLAB Butterworth Low Pass Filter Example 1
% butterworth.m:

% Create a white box on a black background image
M= 256; N = 256;
image = zeros(M,N);
box = ones(64,64);
%box at centre
image(97:160,97:160) = box;

% Show Image

figure
imshow(image)

% compute fft and display its spectra

F=fft2(double(image));
figure
imagesc((abs(fftshift(F))/(M*N)));
colormap(jet);
axis off;
title('FFT Image Spectra');

%compute Butterworth Low Pass Filter

u0 = 20; % set cut off frequency


u=0:(M-1);
v=0:(N-1);
idx=find(u>M/2);
u(idx)=u(idx)-M;
idy=find(v>N/2);
v(idy)=v(idy)-N;
[V,U]=meshgrid(v,u);

for i = 1: M
    for j = 1:N
      %Apply a 2nd order Butterworth  
      UVw = double((U(i,j)*U(i,j) + V(i,j)*V(i,j))/(u0*u0));
      
      H(i,j) = 1/(1 + UVw*UVw);
    end
end    


% display
figure
imshow(fftshift(H))

% Apply filter and do inverse FFT
G=H.*F;
g=real(ifft2(double(G)));

% Show Result

figure
imshow(g)

%
%% MATLAB Low-Pass Butterworth Filter Example 2
% butterworth2.m:

% read in MATLAB demo text image
image = imread('text.png');
[M, N] = size(image);

% Show Image

figure
imshow(image)

% compute fft and display its spectra

F=fft2(double(image));
figure
imshow(abs(fftshift(F))/256);

%compute Ideal Butterworth Pass Filter

u0 = 50; % set cut off frequency


u=0:(M-1);
v=0:(N-1);
idx=find(u>M/2);
u(idx)=u(idx)-M;
idy=find(v>N/2);
v(idy)=v(idy)-N;
[V,U]=meshgrid(v,u);

for i = 1: M
    for j = 1:N
      %Apply a 2nd order Butterworth  
      UVw = double((U(i,j)*U(i,j) + V(i,j)*V(i,j))/(u0*u0));
      
      H(i,j) = 1/(1 + UVw*UVw);
    end
end    
% display
figure
imshow(fftshift(H))

% Apply filter and do inverse FFT
G=H.*F;
g=real(ifft2(double(G)));

% Show Result

figure
imshow(g)

%
%% Low Pass Filtering Noisy Images
%
% How to create noise and results of Low Pass Filtering Use Matlab function,  imnoise() to add noise to image
% See lowpass.m, lowpass2.m examples above.

showDoc = 0;
if showDoc
    doc imnoise;
end

%
%% Other Filters
% Other filters can be created in a similar fashion to above:
% High-Pass Filters— opposite of low-pass, select high frequencies, attenuate those below u_0
% Band-pass— allow frequencies in a range u_0 \ldots u_1 attenuate those outside this range
% Band-reject— opposite of band-pass, attenuate frequencies within u_0 \ldots u_1 select those outside this range
% Notch— attenuate frequencies in a narrow bandwidth around cut-off frequency, u_0
% Resonator— amplify frequencies in a narrow bandwidth around cut-off frequency, u_0
%
%% High Pass Filtering
% A high pass filter is easy to implement and a simple modification on the above Low Pass Filter
% A high Pass Filter is usually defined as 1 - low pass = 1 - H, an ideal High Pass Filter:

% read in MATLAB demo text image
image = imread('text.png');
[M, N] = size(image);


% Show Image

figureheaderdlg;
imshow(image);
title('Original Image');

% compute fft and display its spectra

F=fft2(double(image));
figure;
imagesc((abs(fftshift(F))/(M*N)));
colormap(jet);
axis off;
title('FFT Image Spectra');

%compute Ideal high Pass Filter

u0 = 50; % set cut off frequency


u=0:(M-1);
v=0:(N-1);
idx=find(u>M/2);
u(idx)=u(idx)-M;
idy=find(v>N/2);
v(idy)=v(idy)-N;
[V,U]=meshgrid(v,u);
D=sqrt(U.^2+V.^2);
H=double(D<=u0);

% display
figure;
imshow(fftshift(H));
title('Ideal high Pass Filter');

% Apply filter and do inverse FFT
G=H.*F;
g=real(ifft2(double(G)));

% Show Result

figure;
imshow(g);
title('high Pass Filtered Image');



%%% Demo now to show noisy filtering properties of the high pass filter

% now addnoise to the image

image_noisy = imnoise(double(image),'gaussian');

% show noisy im
figure
imshow(image_noisy)
title('Image with Added Noise');

%%%% compute FFT of noisy image

F = fft2(double(image_noisy));



% Filter using high pass filter H as above

G_highpass=H.*F;

% convert back to image spce with iFFT

g_highpass = real(ifft2(double(G_highpass)));

figure
imshow(g_highpass);
title('high Pass Filtered Noisy Image');



%%%%% Now set Higher cut off frequency and repeat all above

%compute Ideal high Pass Filter

u0 = 80; % set cut off frequency


u=0:(M-1);
v=0:(N-1);
idx=find(u>M/2);
u(idx)=u(idx)-M;
idy=find(v>N/2);
v(idy)=v(idy)-N;
[V,U]=meshgrid(v,u);
D=sqrt(U.^2+V.^2);
H= 1 - double(D<=u0);

% display
figure;
imshow(fftshift(H));
title('High Cut-off Frequency High Pass Filter');

% Apply filter and do inverse FFT


% Filter using High pass filter H as above

G_highpass=H.*F;

% convert back to image spce with iFFT 

g_highpass = real(ifft2(double(G_highpass)));

figure
imshow(g_highpass);
title('High Cut-off Frequency High Pass Filtered Image');

%
%% Convolution
%
% There are many practical, useful applications of convolution – Several important audio and optical effects can be described in terms of convolutions.
% FIltering — In fact the above Fourier filtering is applying convolutions of a low pass filter where the equations are Fourier Transforms of real space equivalents.
% Deblurring — high pass filtering
% Reverb — impulse response convolution (more soon).
% Note: we have seen a discrete real domain example of Convolution with Edge Detection, already.
%
%% 1D Convolution
% Formal Definition of 1D Convolution: Let us examine the concepts using 1D continuous functions.
%
% The convolution of two functions f(x) and g(x), written f(x)*g(x), is defined by the integral:
% f(x)*g(x) = \int_{-\infty}^\infty f(\alpha) g(x-\alpha)\,d\alpha.
% * is the mathematical notation for convolution.
% No Fourier Transform in sight here — but wait!
%
%% 1D Convolution Example in the Real Domain
%
% Convolution of Two Top Hat Functions For example, let us take two top hat functions:
% Let f(\alpha) be the top hat function shown:
% f(\alpha) = 1if |\alpha| \leq 1
%                  0  otherwise,
% and let g(\alpha) be as shown in next slide, defined by:
% f(\alpha) = \frac{1}{2} if 0 \leq |\alpha| \leq 1
%                  0  otherwise,
%
% Our Two Top Hat Functions Plots:
%
% f(\alpha)                                                          g(\alpha)
% The Convolution Process: Graphical Interpretation
% g(-\alpha) is the reflection of this function in the vertical y-axis,
% g(x-\alpha) is the latter shifted to the right by a distance x.
% Thus for a given value of x, f(\alpha)g(x-\alpha) integrated over all \alpha is the area of overlap of these two top hats, as f(\alpha) has unit height.
% An example is shown for x in the range -1\leq x\leq 0 opposite
%
%
% So the solution is: If we now consider x moving from -\infty to +\infty, we can see that
% For x\leq -1 or x\geq 2, there is no overlap;
% As x goes from -1 to 0 the area of overlap steadily increases from 0 to 1/2;
% As x increases from 0 to 1, the overlap area remains at 1/2;
% Finally as x increases from 1 to 2, the overlap area steadily decreases again from 1/2 to 0.
% Thus the convolution of f(x) and g(x), f(x)*g(x), in this case has the form shown on next slide
%
%
% Result of f(x)*g(x)
%
% Mathematically the convolution is expressed by:
% f(x)*g(x) =  (x+1)/2      if -1\leq x\leq 0$},
%                        1/2              if 0\leq x\leq 1, or
%                    1-x/2       if 1\leq x\leq 2$}\\ 0 & \mbox{otherwise.} \end{array} \right.
%
%
%% The Fourier Transform and Convolution
%
% The Convolution Theorem: Convolution in Frequency Space is Easy
% One major reason that Fourier transforms are so important in signal/image processing is the convolution theorem which states that:
%
% If f(x) and g(x) are two functions with Fourier transforms F(u) and G(u), then the Fourier transform of the convolution f(x)*g(x) is simply the product of the Fourier transforms of the two functions, F(u) G(u).
% To obtain the convolution f(x)*g(x)simply compute the inverse Fourier transform of F(u) G(u).
% Recall our Low Pass Filter Example (MATLAB CODE)
%  % Apply filter
%  G=H.*F;
% Formally this is a convolution, where F was the Fourier transform of the image, H the filter
%
%% Computing Convolutions with the Fourier Transform
% Example Applications:
% To apply some reverb to an audio signal.
% To compensate for a less than ideal image capture system.
% More soon --- see Convolution Reverb later in Audio Effects
% Example Application – Deconvolution: Compensating for undesirable imaging effects such a blurring.
% To do this fast convolution we simply:
% Take the Fourier transform of the audio/imperfect image,
% Take the Fourier transform of the function describing the effect of the system,
% Multiply by the effect to apply effect to audio data
% To remove/compensate for effect: Divide by the effect to obtain the Fourier transform of the ideal image.
% Inverse Fourier transform to recover the new improved audio image.
% This process is sometimes referred to as deconvolution.
%
%% Image Deblurring Deconvolution Example
% Image Deblurring Deconvolution Example
% Inverting our Previous Low-Pass Filter Recall our Low Pass (Butterworth) Filter example above: butterworth.m:
% deconv.m and deconv2.m reuses this code and adds a deconvolution stage:
% Our computed butterworth low pass filter,  H is our blurring function
% So to simply invert this we can divide (as opposed to multiply) by  H with the blurred image  G — effectively a high pass filter
%
%                         Ghigh = G./H;
%                         ghigh=real(ifft2(double(Ghigh)));
%                         figure(5)
%                         imshow(ghigh)
%
% In this ideal example we clearly get  F back and to get the image simply to inverse Fourier Transfer.
% In the real world we don't really know the exact blurring function  H so things are not so easy.
%
% deconv.m:

% This used the same script as for butterworth.m example
% Low Pass Filter to create a blurred image
%
% Simply do deconvoution of butterworth low pass (divide rather multiply)
% to invert the result
%
% Note in practice the ideal forward blurring function will not be known 
% SO results will not be so perfect.


% Create a white box on a black background image
M= 256; N = 256;
image = zeros(M,N);
box = ones(64,64);
%box at centre
image(97:160,97:160) = box;

% Show Image

figure
imshow(image)
title('Original Image')

% compute fft and display its spectra

F=fft2(double(image));
figure
imshow(abs(fftshift(F)));
title('Fourier Spectra of Image');

%compute Butterworth Low Pass Filter

u0 = 20; % set cut off frequency


u=0:(M-1);
v=0:(N-1);
idx=find(u>M/2);
u(idx)=u(idx)-M;
idy=find(v>N/2);
v(idy)=v(idy)-N;
[V,U]=meshgrid(v,u);

for i = 1: M
    for j = 1:N
      %Apply a 2nd order Butterworth  
      UVw = double((U(i,j)*U(i,j) + V(i,j)*V(i,j))/(u0*u0));
      
      H(i,j) = 1/(1 + UVw*UVw);
    end
end    


% display
figure
imshow(fftshift(H))
title('FFT Butterworth Low Pass Filter');

% Apply filter and do inverse FFT
G=H.*F;
g=real(ifft2(double(G)));

% Show Result

figure
imshow(g)
title('Butterworth Low Pass Filtered Image');


%%%%%% NOW ILLUSTRATE DECONVOLUTION   %%%%
% Transfer (blurring function) is H --- the low pass filter
% Simply do deconvoution of butterworth low pass (divide rather multiply)
% to invert the result --- effectively now a high pass filter.


Ghigh = G./H;
ghigh=real(ifft2(double(Ghigh)));
figure
imshow(ghigh)
title('Low Pass Image Deconvolved');
 

% High Pass Filter Original Image  Strictly a High Pass Filter is 1 - H 

Fhigh = F.*(1-H);

fhigh=real(ifft2(double(Fhigh)));
figure
imshow(fhigh)
title('Original image High Pass Filtered');

%%% Demo now to show noisy enhancing properties of the High pass filter


% now addnoise to the image

image_noisy = imnoise(image,'gaussian');

% show noisy im
figure
imshow(image_noisy)
title('Image with Noise');


%%%% compute FFT of noisy image

F = fft2(double(image_noisy));



% Filter using high pass filter H as above

G_highpass=F.*(1 - H);

% convert back to image spce with iFFT

g_highpass = real(ifft2(double(G_highpass)));

figure
imshow(g_highpass);
title('High Pass Filter Noisy Image')


% Filter using high pass filter H as above

G_highpass=F./H;

% convert back to image spce with iFFT

g_highpass = real(ifft2(double(G_highpass)));

figure
imshow(g_highpass);
title('Deconvolved Noisy Image')

% Note: Deconvolution is not always as simple as outlined above!  Real noise is harder to deal with and more advance methods (beyond the scope of this module) are needed.
%
% deconv2.m:

% This used the same script as for butterworth2.m example
% Low Pass Filter to create a blurred image
%
% Simply do deconvoution of butterworth low pass (divide rather multiply)
% to invert the result
%
% Note in practice the ideal forward blurring function will not be known 
% SO results will not be so perfect.


% read in MATLAB demo text image
image = imread('text.png');
[M, N] = size(image);

% Show Image

figure
imshow(image)
title('Original Image')

% compute fft and display its spectra

F=fft2(double(image));
figure
imshow(abs(fftshift(F))/256);
title('FFT Spectra of Image');

%compute Ideal Butterworth Pass Filter

u0 = 50; % set cut off frequency


u=0:(M-1);
v=0:(N-1);
idx=find(u>M/2);
u(idx)=u(idx)-M;
idy=find(v>N/2);
v(idy)=v(idy)-N;
[V,U]=meshgrid(v,u);

for i = 1: M
    for j = 1:N
      %Apply a 2nd order Butterworth  
      UVw = double((U(i,j)*U(i,j) + V(i,j)*V(i,j))/(u0*u0));
      
      H(i,j) = 1/(1 + UVw*UVw);
    end
end    
% display
figure
imshow(fftshift(H))
title('Butterworth Low Pass Filter');

% Apply filter and do inverse FFT
G=H.*F;
g=real(ifft2(double(G)));

% Show Result

figure
imshow(g)
title('Butterworth Low Pass Filtered Image');


%%%%%% NOW ILLUSTRATE DECONVOLUTION   %%%%
% Transfer (blurring function) is H --- the low pass filter
% Simply do deconvoution of butterworth low pass (divide rather multiply)
% to invert the result --- effectively now a high pass filter.


Ghigh = G./H;
ghigh=real(ifft2(double(Ghigh)));
figure
imshow(ghigh)
title('Deconvolved Image');

%%% Demo now to show noisy enhancing properties of the High pass filter


% Filter using High pass filter H as above

G_highpass=F./H;

% convert back to image spce with iFFT

g_highpass = real(ifft2(double(G_highpass)));

figure
imshow(g_highpass);
title('Origial Image Deconvolved');



% High Pass Filter Original Image  Strictly a High Pass Filter is 1 - H 

Fhigh = F.*(1-H);

fhigh=real(ifft2(double(Fhigh)));
figure
imshow(fhigh)
title('Original image High Pass Filtered');

% now addnoise to the image

image_noisy = imnoise(double(image),'gaussian');

% show noisy im
figure
imshow(image_noisy)
title('Image with Noise Added');


%%%% compute FFT of noisy image

F = fft2(double(image_noisy));

% Filter using high pass filter H as above

G_highpass=F.*(1 - H);

% convert back to image spce with iFFT

g_highpass = real(ifft2(double(G_highpass)));

figure
imshow(g_highpass);
title('High Pass Filter Noisy Image');


% Filter using high pass filter H as above

G_highpass=F./H;

% convert back to image spce with iFFT

g_highpass = real(ifft2(double(G_highpass)));

figure
imshow(g_highpass);
title('Deconvolved Noisy Image');




% now addnoise to the blurred low pass image g and do the same

image_noisy = imnoise(double(g),'gaussian');

% show noisy im
figure
imshow(image_noisy)
title('Low Pass Filtered (Blurred) Image with Noise')


%%%% compute FFT of noisy image

F = fft2(double(image_noisy));

% Filter using deconvolution H as above

G_highpass=F./H;

% convert back to image spce with iFFT

g_highpass = real(ifft2(double(G_highpass)));

figure
imshow(g_highpass);
title('Deconvolved Blurred Noisy Image');



% Filter using High Pass Filter (1 - H) as above

G_highpass=F.*(1 - H);

% convert back to image spce with iFFT

g_highpass = real(ifft2(double(G_highpass)));

figure
imshow(g_highpass);
title('High Pass Filtered Blurred Noisy Image');

%
%% A Practical Audio Example of Convolution: Reverb
%
% Reverberation Reverb (for short) is probably one of the most heavily used effects in audio.
% Reverberation is the result of the many reflections of a sound that occur in a room.
% From any sound source, say a speaker of your stereo, there is a direct path that the sounds covers to reach our ears.
%         The Spaciousness of a Room:
% Sound waves can also take a slightly longer path by reflecting off a wall or the ceiling, before arriving at your ears.
%
% A reflected sound wave like this will arrive a little later than the direct sound, since it travels a longer distance, and is generally a little weaker, as the walls and other surfaces in the room will absorb some of the sound energy.
% Reflected waves can again bounce off another wall before arriving at your ears, and so on.
% This series of delayed and attenuated sound waves is what we call reverb, and this is what creates the spaciousnesssound of a room.
% Clearly large rooms such as concert halls/cathedrals will have a much more spaciousness reverb than a living room or bathroom.
%% Reverb Simulation
% How to implement a Digital Audio Effect There are a few ways to simulate reverb.
%
% There are many ways to simulate reverb.
% Two Broad classes of approach studied here (there are others):
% Filter Bank/Delay Line methods
% Convolution/Impulse Response methods
%% Schroeder's Reverberator
% Schroeder's model of Reverb (1961)
% Early digital reverberation algorithms tried to mimic the a rooms reverberation by using primarily consisted of two types of infinite impulse response (IIR) filters with the aim to make the output  gradually decay.
%
% Historical note: Much of the early work on digital reverb was performed by Schroeder.
%
%% Schroeder's reverberator example (one of a few variations):
%
% This particular design uses four comb filters and two allpass filters:
%
% Note: This design does not create the increasing arrival rate of reflections, and is rather primitive when compared to current algorithms.
%
%% A simple Schroeder reverberator,  schroeder1.m:
% We basiclly connect n allpass filters in series.
% No Comb Filters yet.
% function [y,b,a]=schroeder1(x,n,g,d,k)
% %This is a reverberator based on Schroeder's design which consists of n all pass filters in series.
% %
% %The structure is:  [y,b,a] = schroeder1(x,n,g,d,k)
% %
% %where x = the input signal
% %      n = the number of allpass filters
% %      g = the gain of the allpass filters (this should be less than 1 for stability)
% %      d = a vector which contains the delay length of each allpass filter
% %      k = the gain factor of the direct signal
% %      y = the output signal
% %      b = the numerator coefficients of the transfer function
% %      a = the denominator coefficients of the transfer function
% %
% % note: Make sure that d is the same length as n.
%
%
% % send the input signal through the first allpass filter
% [y,b,a] = allpass(x,g,d(1));
%
% % send the output of each allpass filter to the input of the next allpass filter
% for i = 2:n,
%    [y,b1,a1] = allpass(y,g,d(i));
%    [b,a] = seriescoefficients(b1,a1,b,a);
% end
%
% % add the scaled direct signal
% y = y + k*x;
%
% % normalize the output signal
% y = y/max(y);
%
% end
%
% The support files to do the filtering (for following reverb methods also) are linked here and listed in the supporting functions below :
% delay.m,
% seriescoefficients.m,
% allpass.m
%
% Example call, reverb_schroeder_eg.m:

% read the sample waveform
filename = 'acoustic.wav';
if exist(filename) == 0
    % read file from URL
    [monox,Fs] = webread('http://users.cs.cf.ac.uk/Dave.Marshall/Multimedia/Lecture_Examples/Digital_Audio_FX/acoustic.wav');
else
[monox,Fs] = audioread(filename);
end

% Call Schroeder1 reverb
%set the number of allpass filters
n = 6;
%set the gain of the allpass filters
g = 0.9;
%set delay of each allpass filter in number of samples
%Compute a random set of milliseconds and use sample rate

s = rng;
d = floor(0.05*rand([1,n])*Fs);
%set gain of direct signal
k= 0.2;


[y, b, a] = schroeder1(monox,n,g,d,k);

% The input signal (blue) and reverberated signal (red):

figure;
hold on
plot(y,'r');
plot(monox,'b');
title('Schroeder Reverberated Signal');

%

% Play schroeder1 Reverb Audio
sound_choice = y;

sound(sound_choice,Fs);

 

% write output
%audiowrite('out_schroederreverb.wav', y,Fs);

%
%
%% Convolution Reverb Implementation
% Commercial Convolution Reverbs
% Altiverb — one of the first mainstream convolution reverb effects units
%
% Most sample based synthesisers (E.g. Kontakt, Intakt) provide some convolution reverb effect
% Dedicated sample-based software instruments such as Garritan Violin and PianoTeq Piano use convolution not only for reverb simulation but also to simulate key responses of the instruments body vibration.
%
% All major digital audio workstaions provide a convolution reverb plugin.
%
% These days a lot of audio software include some convolution reverb facility.
%% Room Impulse Responses
% Record a Room Impulse Apart from providing a high (professional) quality recording of a room's impulse response, the process of using an impulse response is quite straightforward:
% Record a short impulse (gun shot,drum hit, hand clap) in the room.
% Room impulse responses can be simulated in software also.
% The impulse encodes the rooms reverb characteristics:
%
%% MATLAB Convolution Reverb
% Let's develop a fast convolution routine: fconv.m
% function [y]=fconv(x, h)
% %FCONV Fast Convolution
% %   [y] = FCONV(x, h) convolves x and h, and normalizes the output
% %         to +-1.
% %
% %      x = input vector
% %      h = input vector
% %
% %
%
% Ly=length(x)+length(h)-1;  %
% Ly2=pow2(nextpow2(Ly));    % Find smallest power of 2 that is > Ly
% X=fft(x, Ly2);		   % Fast Fourier transform
% H=fft(h, Ly2);	           % Fast Fourier transform
% Y=X.*H;        	           %
% y=real(ifft(Y, Ly2));      % Inverse fast Fourier transform
% y=y(1:1:Ly);               % Take just the first N elements
% y=y/max(abs(y));           % Normalize the output
%
% end
%
% See also: MATLAB built in function conv()

help conv;

%
% A MATLAB example of convolution reverb in action, reverb_convolution_eg.m:

% read the sample waveform
filename = 'acoustic.wav';
if exist(filename) == 0
    % read file from URL
    [x,Fs] = webread('http://users.cs.cf.ac.uk/Dave.Marshall/Multimedia/Lecture_Examples/Digital_Audio_FX/acoustic.wav');
else
[x,Fs] = audioread(filename);
end

% read the impulse response waveform

% read the sample waveform
filename = 'impulse_demo.wav';
if exist(filename) == 0
    % Download Impulse_Responses
    unzip('http://users.cs.cf.ac.uk/Dave.Marshall/Multimedia/Lecture_Examples/Digital_Audio_FX/Reverb/Impulse_Responses.zip');
    addpath('Impulse_Responses/')

end

Impulse_filename= 'impulse_MM.wav';


[imp,Fsimp] = audioread(Impulse_filename);


% Do convolution with FFT
y = fconv(x,imp);


figure;
hold on
plot(x,'b');
plot(y,'r');
title('Impulse Response Reverberated Signal');


figure;

plot(imp);
title('Impulse Response');

%

% Play Convolution Reverb Audio also able to select/play the impulse
% response
sound_choice = y;

sound(sound_choice,Fs);

 

% write output
%audiowrite('out_IRreverb.wav', y,Fs);

% Note: It is easy to implement some other (odd non-reverb) effects also:
% Select the "Reverse Cathedral Impulse Response" Convolution Reverb above.
% Select the "Speech Impulse Response Convolution Reverb!"
%
%% Further Reading
% DAFX: Digital Audio Effects, Udo Zolzer, John Wiley and Sons Ltd , 2002 (ISBN-13: 978-0471490784)
%  Excellent coverage of audio, signal processing effects and synthesis plus a lot more
%  All MATLAB examples
%  ONLINE Copies in library
%
% Digital Image Processing Using MATLAB, Rafael C. Gonzalez, Richard E. Woods and Steven L. Eddins, Prentice Hall, 2004  (ISBN-13: 978-0130085191)
%
% Excellent coverage of Image processing examples plus a lot more
% All MATLAB examples
% Useful for CM0311 Image Processing
% Copies in library
%
% Additional Notes/Examples:
% Audio Synthesis— CM3106 Multimedia Module Notes
% Audio Effects— CM3106 Module Notes
%
%% Supporting functions
%
% schroeder1.m:

function [y,b,a]=schroeder1(x,n,g,d,k)
%This is a reverberator based on Schroeder's design which consists of n all pass filters in series.
%
%The structure is:  [y,b,a] = schroeder1(x,n,g,d,k) 
%
%where x = the input signal
%      n = the number of allpass filters	
%      g = the gain of the allpass filters (this should be less than 1 for stability)
%      d = a vector which contains the delay length of each allpass filter
%      k = the gain factor of the direct signal
%      y = the output signal
%      b = the numerator coefficients of the transfer function
%      a = the denominator coefficients of the transfer function
%
% note: Make sure that d is the same length as n.


% send the input signal through the first allpass filter
[y,b,a] = allpass(x,g,d(1));

% send the output of each allpass filter to the input of the next allpass filter
for i = 2:n
   [y,b1,a1] = allpass(y,g,d(i));
   [b,a] = seriescoefficients(b1,a1,b,a);
end   

% add the scaled direct signal
y = y + k*x;

% normalize the output signal
y = y/max(y);

end

%
% delay.m:

function y=delay(x,d)

%This function delays a signal by d samples.
%
%The structure is: y = delay(x,d)
%
%where x = input signal
%      d = the number of samples to delay the input signal
%


y = [zeros(d,1); x];

end

% seriescoefficients.m:

function [b,a]=seriescoefficients(b1,a1,b2,a2)

%This function gives the filter coefficients of the series connection of two filters.
%
%The structure is:  [b,a] = seriescoefficients(b1,a1,b2,a2)
%
%where b1 = the numerator coefficients of the 1st transfer function
%      a1 = the denominator coefficients of the 1st transfer function
%      b2 = the numerator coefficients of the 2nd transfer function
%      a2 = the denominator coefficients of the 2nd transfer function
%      b = the numerator coefficients of the composite transfer function
%      a = the denominator coefficients of the composite transfer function
%

b = conv(b1,b2);

a = conv(a1,a2);

end

%
% allpass.m:

function [y,b,a]=allpass(x,g,d)

%This is an allpass filter function.
%
%The structure is:  [y,b,a] = allpass(x,g,d)
%
%where x = the input signal
%      g = the feedforward gain (the feedback gain is the negative of this) (this should be less than 1 for stability)
%      d = the delay length
%      y = the output signal
%      b = the numerator coefficients of the transfer function
%      a = the denominator coefficients of the transfer function
%



%If the feedback gain is more than 1, set it to 0.7 .
if g>=1
   g=0.7;
end   

%Set the b and a coefficients of the transfer function depending on g and d.
b=[g zeros(1,d-1) 1];
a=[1 zeros(1,d-1) g];

%filter the input signal 
y=filter(b, a,x);

end

%
%  fconv.m

function [y]=fconv(x, h)
%FCONV Fast Convolution
%   [y] = FCONV(x, h) convolves x and h, and normalizes the output  
%         to +-1.
%
%      x = input vector
%      h = input vector
% 
%  

Ly=length(x)+length(h)-1;  % 
Ly2=pow2(nextpow2(Ly));    % Find smallest power of 2 that is > Ly
X=fft(x, Ly2);		   % Fast Fourier transform
H=fft(h, Ly2);	           % Fast Fourier transform
Y=X.*H;        	           % 
y=real(ifft(Y, Ly2));      % Inverse fast Fourier transform
y=y(1:1:Ly);               % Take just the first N elements
y=y/max(abs(y));           % Normalize the output

end