%=========================================================
%  Program  : LS_DFT.m
%  Author   : L. Sevgi
%  Date     : May 2006
%  Purpose  : To take the FFT of two sinuzoids
%=========================================================
clear all;
close all;
clc;
%  Get the input parameters
     fr1=input('Frequency of the first sinusoid [Hz] = ? ');
     a1=1; %input('Amplitude of the first sinusoid [ ] = ? ');
     fr2=input('Frequency of the second sinusoid [Hz] = ? ');
     a2=0.5; %input('Amplitude of the first sinusoid [ ] = ? ');
     T=input ('Time record length  [s] = ? ');
     dt=input ('sampling time interval  [s] = ? ');
     fmax=input('maximum frequency for DFT [Hz] = ? ');
     df=input('frequency sampling interval for DFT [Hz] = ? ');

     N=T/dt;
     w1=2*pi*fr1;
     w2=2*pi*fr2;
     M=fmax/df;
    
% build the input time series
    for k=1:N
       st(k)=a1*sin(w1*dt*k)+a2*sin(w2*dt*k);
%      st(k)=a1*exp(j*w1*dt*k)+a2*exp(j*w2*dt*k);
    end
% Apply the DFT with M points
    for k=1:M 
        Sf(k)=complex(0,0);   
        for n=1:N
            Sf(k)=Sf(k)+st(n)*exp(-i*2*pi*n*dt*k*df);
        end
        Sf(k)=Sf(k)*dt;
    end
% Prepare the frequency samples
    for k=1:M
        F(k)=(k-1)*df;
    end
% plot the output
    plot(F,abs(Sf));
    title('The DFT of the Sum of two Sinusoids')
    xlabel('Frequency [Hz]')
    ylabel('Amplitude')
% --------------- END ---------------