% 1-D DFT example
% fft algorithm and formula analysis
clc; close all; clear all;

% take a complex sinusoidal pattern for example
% f(t)=10+15*sin(2*pi*50*t+pi/3)-20*cos(2*pi*75-pi/6);

%%  fft algorithm
% component frequencies in the sigal, known in prior
% f0=0; dc component
f1=50;
f2=80;

% define sampling interval (detT) and time span
Fs=250; % frequency of sampling, detT=1/Fs;
detT=1/Fs;
Ns=200; % number of sampling, time duration of sampling T=(Ns-1)*detT
T=(Ns-1)*detT;
tspan=(0:Ns-1)*detT;
fval=10+15*sin(2*pi*f1*tspan+pi/3)-20*cos(2*pi*f2*tspan-pi/6); % discrete values

% define the fequency numbers and interval (detF) to be computed by fft
Nf=Ns; % equal to the number of sampling points
% Nf=240; % other values, may cause aliasing frequency in the frequency spectrum
detF=1/T; % detF=1/T=Fs/(Nf-1), so the frequency range is [0, Fs]
fspan=(0:Nf-1)*detF; % 0:Fs/(Nf-1):Fs
% The maximum frequency is  the sampling frequency Fs

% fft computation
Fcoef=fft(fval,Nf); % Fcoef: complex, Fourier transform coefficients
Fcoef=abs(Fcoef); % Fourier spectrum magnitude
% magnitude adjustment
Fmag=Fcoef/(Nf/2); % two symmetric parts are added together
Fmag(1)=Fcoef(1)/Nf; % the be beginning value alone for dc

% display;
figure;subplot(211);
plot(tspan,fval,'r-','linewidth',1.5); title('raw pattern');
hold on; plot(tspan,fval,'b*','MarkerSize',6);
subplot(212);
% only the first half part of the frequency spectrum is displayed
stem(fspan(1:round(Nf/2)),Fmag(1:round(Nf/2))); title('Frequency spectrum');

% Comments: the more number Ns, the more chance to retrieve component frequency
% when Nf~=Ns,the positions of component frequency will be inaccurate

%%　define a function for direction 1-D dft

function=1D_dft(f,u)

N=length(f);

if nargin<2

u=0:N-1

end

F=0

omega=2*pi*u/N; % angular frequency

for  x=0:N-1 % spatial or time postions

F=F+f(x+1)*cos(omega*x)-1i*sin(omega*x)); % cumulative sum

end

F=F/N; % normalization

return

%% application of DFT formula
% Given, N sampling points that are regularly spaced over the duration
% Fcoef(u)=(1/N)*sum(fval*exp(-2i*pi*u*x/N))
% where x=0,1,2,...,N-1 and u=0,1,2,...,N-1
% or, Fcoef(u)=(1/N)*sum(fval*W^(-uk)), W=exp(2*pi*i/N)

N=numel(fval);
x=0:N-1; % x-domain coordinates, f(x)=fval, not tspan
W=exp(2*1i*pi/N);
Fcoef=ones(1,N);

% define frequency sequencies
uspan=linspace(0,Fs,N); % but in real cases, Fs is unknown in general

for u=0:N-1; % frequency sequence
   Fcoef(u+1)=sum(fval.*W.^(u*x))/N; % divided by N
end
Fcoef=abs(Fcoef);
Fmag=2*Fcoef;
Fmag(1)=Fcoef(1);
% magnitude ajustment

u=0:N-1;

figure; subplot(211); plot(x,fval); title('raw pattern');
subplot(212); stem(uspan(1:round(N/2)),Fmag(1:round(N/2))); title('Frequency spectrum');
% Comments: the position of component frequency is inaccurate
% if the frequency interval are not defined accurately