% Filename: FrequencyDomainFiltering.m
% Author: Andreas Schwarzinger                                                Date: 11/24/2019
% Source: www.dsp-vault.com
%
% Notes: 
% This script illustrates frequency domain filtering via two separate methods.
% 1. The first method splits the input waveform, x[n], into non-overlapping sections,
%    takes their DFT, multiplies the results by H[m], and takes the IDFT of each product
%    to arrive at the non-overlapping output sections that are then concatenated to
%    form the final output signal y[n].
% 2. The second method is very similar but divides up the input waveform x[n] into 
%    overlapping section, where the start and end of each section is terminated by
%    a special mask that forces the first and last sample toward zero.

clc;                % clear the command window
clear variables;    % clear variables from a previous run
close all;          % close all figures

% ------------------------------------------------------------------------------
% 0. Simulation Setup
Fs  = 1e6;                 % Sample Rate in Hz
Len = 4*4096;              % Number of samples in test waveform 
n   = 0:Len-1;             % The sample indices
disp(['Nyquist range spans: ' num2str(-Fs/2) ' to ' num2str(Fs/2) ' Hz.']);
 
 
% -------------------------------------------------------------------------------
% 1a. Generate input signal 1 -> cos(2*pi*n*F1/Fs)
F1      = 3e3;                  % Frequency of test signal 1
Signal1 = cos(2*pi*n*F1/Fs);    % Signal1
 
% 1b. Generate input signal 2 -> cos(2*pi*n*F2/Fs)
F2      = 10e3;                 % Frequency of test signal 2
Signal2 = cos(2*pi*n*F2/Fs);    % Signal2

% 1c. Generate input signal 3 -> cos(2*pi*n*F3/Fs)
F3      = 75e3;                 % Frequency of test signal 3
Signal3 = cos(2*pi*n*F3/Fs);    % Signal3
 
% 1d. The final composite signal x[n]
x_n = Signal1 + Signal2 + Signal3;
 
% In this exercise we will be filtering away Signal3. To evaluate how well
% the filtering process worked, we synthesize the ideal filtered signal.
x_n_Filtered = Signal1 + Signal2;
 
% Let's take a look at the composite waveform, x[n]
figure(1);
DisplayRange = 3900:4200;
subplot(3,1,1);
plot(DisplayRange, x_n(1, DisplayRange), 'k'); grid on; hold on;
plot(DisplayRange, x_n_Filtered(1, DisplayRange), 'k.-', 'LineWidth', 2);
title('Input Sequences');
legend('Input Waveform, x[n]', 'Ideal Filtered Waveform');


% ----------------------------------------------------------------------------
% 2. Defining the Filter frequency Response. Note, the frequency response will
% be entirely real valued.
% Our goal is the reject the signal with Frequency F3.
% Remember, that the frequency step is equal to Fs/SectionLength
NumSections        = 4;
SectionLength      = 4096;
FStep              = Fs/SectionLength;
disp(['Frequency Step: ' num2str(FStep)]);

H_m                = zeros(1, SectionLength);
H_m(1,1:160)        = ones(1, 160);   % Define passband for positive frequencies
H_m(1,end-158:end)  = ones(1, 159);   % Define passband for the negative frequencies

% ----------------------------------------------------
% 3a. Method 1: Non-Overlapping sections
 
y1_n = zeros(1, length(x_n));
for SectionIndex  = 0:(NumSections-1)
   Range          = (1:SectionLength) + SectionIndex*SectionLength; 
   Section        = x_n(1, Range);
   DFT_Section    = fft(Section);
   Y_m            = H_m .* DFT_Section;
   y1_n(1, Range) = ifft(Y_m);
   if(SectionIndex == 0)
     figure(2);
     m = 0:600;
     f = m*FStep;
     % The division by 500 is only there so that we can easily compare
     % the signal spectrum and the filter's magnitude response.
     plot(f, abs(DFT_Section(1, m+1))/500, 'k'); grid on; hold on; 
     plot(f, abs(H_m(1, m+1)), 'k-', 'LineWidth', 2);
     title('Signal Spectrum and Filter Response');
     xlabel('Frequency in Hz');
     legend('Signal Spectrum', 'Magnitude Response of Filter');
   end
end

% ---------------------------------------------------
% 3b. Method 2: Overlapping Sections
 
HannLength    =  256;               % 256 samples long – try 512 or 128 if you like
HanningWindow = hann(HannLength)';
FrontMask  = HanningWindow(1, 1:HannLength/2);               % 128 samples long
BackMask   = HanningWindow(1, (1+HannLength/2):HannLength);  % 128 samples long
TotalMask  = [FrontMask, ones(1,SectionLength-HannLength), BackMask]; % 4096 samples
y2_n       = zeros(1, length(x_n));
for SectionIndex  = 0:(NumSections-1)
   Range          = (1:SectionLength) + SectionIndex*(SectionLength - HannLength/2); 
   Section        = x_n(1, Range).*TotalMask;
   DFT_Section    = fft(Section);
   Y_m            = DFT_Section .* H_m ;
   OutputSection  = ifft(Y_m);
   y2_n(1, Range) = y2_n(1, Range) + OutputSection;
end
 
figure(1);
subplot(3,1,2);
plot(DisplayRange, y1_n(1, DisplayRange), 'k.-');  grid on; hold on;
plot(DisplayRange, y2_n(1, DisplayRange), 'k');
axis([3900 4200 -2 3]);
title('Filtered Waveforms via Method1 and Method2');
legend('Method1', 'Method2');

subplot(3,1,3);
plot(DisplayRange, y2_n(1, DisplayRange) - x_n_Filtered(1, DisplayRange), 'k'); grid on; hold on;
%plot(DisplayRange, y2_n(1, DisplayRange) - x_n_Filtered(1, DisplayRange), 'k'); grid on; hold on;
title('Error between Ideally Filtered Fignal and Method2 Result.');
 
 
figure(3);
subplot(2,1,1);
stem(FrontMask, 'k'); grid on;
title('The front part of the Hanning Window Magnitude Mask');
subplot(2,1,2);
stem(BackMask, 'k'); grid on;
title('The back part of the Hanning Window Magnitude Mask');
