%============================================================
% Module : dft_ex.m
% Description : Test if Dtft algorithm works well.
% History : 2002/10/13
%           Brian(CheolHoon) Park
%
% for m=1:SampleNo;
%	w=2*pi*m*HarmHz/Fs;
%	z=exp(-j*w); 		% z = cos(w) - j*sin(w) = cos(2*pi*m/SecNum) -j*sin(2*pi*m/SecNum); 
%	
%	X=0;
%	for n = 1:SampleNo
%		X = x(n+1) + z. * X; %Xr2=  x(n) + cos(2*pi*m/SecNum) *Xr + sin(2*pi*m/SecNum) *Xi;
%							 %Xi2=         cos(2*pi*m/SecNum) *Xi - sin(2*pi*m/SecNum) *Xr;
%  	end;
% 	MagX(m)=abs(X)*(2/SampleNo);
% end;
%=============================================================
clear all;
close all;

HarmHz=120;
SecNum=232;
Fs=HarmHz*SecNum;
Ts=1/Fs;
REV_N=1;
SampleNo=SecNum*REV_N;
t=0:Ts:(SampleNo-1)*Ts;
SecX=1:SampleNo;
n=0:SampleNo-1;
ADC_Scale=2^8;

fn=SecNum;			 %fs/harmoincs (x times of HarmHz)
fs=SecNum*HarmHz;		 %sampling freq(Hz)
freq_t	= fn*(0:SampleNo/2-1)/SampleNo;    %if plot with harmoincs
%freq_t	= fs*(0:SampleNo/2-1)/SampleNo;    %if plot with Hz

f= 10*HarmHz;
%x=round(10*sin(2*pi*f*t));
%x=10*sin(2*pi*f*t);
x=10*sin(2*pi*f*t)+5*cos(pi*f*t);

%====================================================================
% Hamming window
%====================================================================
w=2*(0.54-0.46*cos(2*pi*n/(SampleNo-1)));     
xL=w.*x;

%====================================================================
% DTFT of x
%====================================================================
for m=1:SampleNo;
       CosTbl(m)= round(ADC_Scale*cos(2*pi*(m-1)/SecNum));
       SinTbl(m)= round(ADC_Scale*sin(2*pi*(m-1)/SecNum));
end;

for m=1:SampleNo;

	Xr=0;
	Xi=0;
	for n = SampleNo:-1:1
	   	Xr2=  xL(n) + round(round((CosTbl(m)*Xr + SinTbl(m)*Xi))/ADC_Scale);
	   	Xi2=          round(round((CosTbl(m)*Xi - SinTbl(m)*Xr))/ADC_Scale);
		Xr=Xr2;
		Xi=Xi2;
	end;
	MagX(m)=sqrt(Xr^2 + Xi^2)*(2/SampleNo);
	%MagX(m)= (Xr^2 + Xi^2)*(4/SampleNo^2);
end;

%====================================================================
% DTFT of xL
%====================================================================
for m=1:SampleNo;
	CosTbl(m)= round(ADC_Scale*cos(2*pi*(m-1)/SecNum));
	SinTbl(m)= round(ADC_Scale*sin(2*pi*(m-1)/SecNum));
end;

for m=1:SampleNo;

	Xr=0;
	Xi=0;
	for n = SampleNo:-1:1
	   	Xr2=  x(n) + round(round((CosTbl(m)*Xr + SinTbl(m)*Xi))/ADC_Scale);
	   	Xi2=         round(round((CosTbl(m)*Xi - SinTbl(m)*Xr))/ADC_Scale);
		Xr=Xr2;
		Xi=Xi2;
	end;
	MagXL(m)=sqrt(Xr^2 + Xi^2)*(2/SampleNo);
	%MagXL(m)= (Xr^2 + Xi^2)*(4/SampleNo^2);
end;


%====================================================================
% FFT
%====================================================================
scale_t=2/SampleNo;
xm=x-mean(x);	     %remove DC	
x_fft=scale_t*fft(xm);
X_FFT=abs(x_fft);
[max_mag,max_i]=max(X_FFT(1:SampleNo/2));

%====================================================================
% DFT of x
%====================================================================
for iFreq=0:SampleNo-1
    x_sin=round(ADC_Scale*sin(2*pi*iFreq*HarmHz*t));
    x_cos=round(ADC_Scale*cos(2*pi*iFreq*HarmHz*t));
    sum_sin=round(sum(x.*x_sin)/ADC_Scale);
    sum_cos=round(sum(x.*x_cos)/ADC_Scale);
    c_sin(iFreq+1)=2*sum_sin/(SampleNo);
    c_cos(iFreq+1)=2*sum_cos/(SampleNo);
    X_DFT(iFreq+1)=sqrt(c_sin(iFreq+1)^2 + c_cos(iFreq+1)^2);
end;

%====================================================================
% DFT of xL
%====================================================================
for iFreq=0:SampleNo-1
    x_sin=round(ADC_Scale*sin(2*pi*iFreq*HarmHz*t));
    x_cos=round(ADC_Scale*cos(2*pi*iFreq*HarmHz*t));
    sum_sin=round(sum(xL.*x_sin)/ADC_Scale);
    sum_cos=round(sum(xL.*x_cos)/ADC_Scale);
    c_sin(iFreq+1)=2*sum_sin/(SampleNo);
    c_cos(iFreq+1)=2*sum_cos/(SampleNo);
    XL_DFT(iFreq+1)=sqrt(c_sin(iFreq+1)^2 + c_cos(iFreq+1)^2);
end;


x_axis=freq_t(1:SecNum/2);
xHarmHz = 0:SampleNo-1;

figure(1);
	plot(xHarmHz,x,xHarmHz,w);
	zoom on; axis on; grid on;

figure(2);
	plot(1:116,MagX(1:length(x_axis)),'b',1:116,MagXL(1:length(x_axis)),'r');
	zoom on; axis on; grid on; legend('w/o window','w/ window');

figure(3);
	plot(x_axis,X_FFT(1:length(x_axis)));
	zoom on; axis on; grid on;

figure(4);
	plot(1:116,X_DFT(1:length(x_axis)),'b',1:116,XL_DFT(1:length(x_axis)),'r');
	%plot(1:116,X_DFT(1:length(x_axis)));
	zoom on; axis on; grid on; legend('w/o window','w/ window');