%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% This is a code to show how a FT match filters out signals at a given frequency. 
% All that a Fourier Transform does is compute the average of the time-domain 
% signals after having rotated the phases of the complex vector at each point in 
% time appropriately. The amount of rotation is "Period-dependent". If P is the 
% period, then the rotation at each time (assuming t=0 as the reference) would be 
% 
%                                                -- wr, 21 Sep 2017 
% ========================================================
% Code is written to work in Octave. 
% You need to have teh geometry package from octave-forge 
% for the arrows to show. 
% To get the geometry package, from the octave prompt do 
% the following: 
%       >pkg install -forge geometry 
%       >pkg load geometry 
% 
% You also need GhostScript for saving the figures. 
% ========================================================
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

clc 
clear all
%close all

N = 64;   % Number of time-samples 
f1 = 0.3; % Frequency of Input signal 1
f2 = 1;   % Frequency of Input signal 2
A1 = 1.;  % Amplitude of Input signal 1 
A2 = 0.0; % Amplitude of Input signal 2 

t = linspace(0,2*pi,N);   % Time array 
s1 = A1*exp(i*2*pi*f1*t); % Signal at frequency f1 
s2 = A2*exp(i*2*pi*f2*t); % Signal at frequency f2 
s_in = s1 + s2;           % A simple linear combination of 

% Pick a reference frequency at which you wish to 
% measure the power: 
f_ref = f2;               % 
s_ref = exp(i*2*pi*f_ref*t) ; % The reference signal to find the FT at. 
                              % Use a signal with monochromatic frequency. 
			      % Here I use either f1 or f2. 

% =======================================================
A = mean(abs(s_in)); 

%++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% How does the signal look like in the Time domain: 
subplot(2,1,1),plot(t,real(s_in),'r-o',t,imag(s_in),'g+') 
legend('real','imag')
%legend('real','imag','location','north')

%++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% Some careful setting to make the arrows discernible: 
scale = 0.3; % A factor to control the length of the arrows. 
             % use this together with the ratio "r". 
r = A*scale;
L = A*ones(1,N);
W = 0.1*ones(1,N);
R = 0.1*ones(1,N);
H = 0.5*ones(1,N); 
% Get the tail and head coord of the arrows: 
x1 = t; 
y1 = abs(s_in); 
th1 = arg(s_in); % Angle 
x2 = x1 + r*cos(th1); 
y2 = y1 + r*sin(th1); 
orig = [[x1;y1]'];
tip = [[x2;y2]'];
arrows = [orig tip];
drawArrow(arrows, L, W, R, H)
axis([0,2*pi,-2.2,2.2])
xlabel('TIME')
ylabel ('Amplitude')
title_str = strcat('Input Signal at f1: ',num2str(f1),'Hz, A1: ',num2str(A1), '; f2: ',num2str(f2), ',Hz, A2: ',num2str(A2));
title(title_str)

% Rotate the phases (Demonstrating the DFT using simple geometry);  
th = arg(s_ref); % The reference signal you wish to de-rotate your phases by
for j = 1:N 
	th_now = -th(j) ; 
	Rota = [cos(th_now) -sin(th_now); sin(th_now) cos(th_now)];
	re = real(s_in(j)); 
	im =imag(s_in(j)); 
	C = [re; im]; 
	C = Rota*C;
	re = C(1);
	im = C(2); 
	s_out(j) = re + i*im; 
end
th = arg(s_out)+pi/2; 
x1 = t; 
y1 = abs(s_out); 
A = mean(y1); % This is what the output power would be at the reference frequency. 
r = A*scale;

L = A*ones(1,N);
W = 0.1*ones(1,N);
R = 0.1*ones(1,N);
H = 0.5*ones(1,N); 
x2 = x1 + r*cos(th); 
y2 = y1 + r*sin(th); 
orig = [[x1;y1]'];
tip = [[x2;y2]'];
arrows = [orig tip];
subplot(2,1,2),plot(t,real(s_out),'r-o',t,imag(s_out),'g+') 
drawArrow(arrows, L, W, R, H)
axis([0,2*pi,-2.2,2.2])
power = mean(real(s_out))
xlabel('TIME')
ylabel ('Amplitude')
title_str = strcat('Output Signal extracted at f: ',num2str(f_ref),'Hz, A: ',num2str(power));
title(title_str)
figure(1) 
outfile = strcat('Input_',num2str(f1),'Hz_',num2str(A1),'_',num2str(f2),'Hz_',num2str(A2),'_Output_at_',num2str(f_ref),'Hz.png')

%print -djpg outfile.jpg 
print(outfile,'-dpng')
