I read this DFT code some where,
Get the input parameters
fr1=input('Frequency of the first sinusoid [Hz] = ? ');
a1=input('Amplitude of the first sinusoid [ ] = ? ');
fr2=input('Frequency of the second sinusoid [Hz] = ? ');
a2=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); % Here they are calculating time series
end
% Apply the DFT with M points
for k=1:M
Sf(k)=complex(0,0); % I didnt get why this step
for n=1:N
Sf(k)=Sf(k)+st(n)*exp(i*2*pi*n*dt*k*df);
end
Sf(k)=Sf(k)*dt/2*pi; % And why this one
end
% Prepare the frequency samples
for k=1:M
F(k)=(k1)*df;
end
% Plot the output
plot(F,abs(Sf));
title('The DFT of the Sum of two Sinusoids')
xlabel('Frequency [Hz]'); ylabel('Amplitude')
This is from Levent Sevgi
DOGUS University, Electronics and Communication Engineering Department
Zeamet Sokak, No 21, Acıbadem, Istanbul, Turkey, lsevgi@dogus.edu.tr
"Wayne King" <wmkingty@gmail.com> wrote in message <hlolgh$eca$1@fred.mathworks.com>...
> "kk KKsingh" <akikumar1983@gmail.com> wrote in message <hlo4ke$4on$1@fred.mathworks.com>...
> > I want to write my own code for DFT for regular sampling signall
> >
> > Expression looks like this
> >
> > Objective  Apply DFT to go from X domain to wave number domain K! which is simlar to go from t domain to frequency domain
> >
> >
> > %P(k,omega)=deltax sum{n=0,N1} P(ndeltax,omega)exp(j k n deltax)
> >
> > WHERE
> > deltax= sampling interval
> > n=number of samples we over which we are doing summation
> > M=number of frequency domain sample assuming N=M
> > k=m*deltak
> > m=M/2:M/21
> > deltak=2*pi/deltak
> >
> >
> >
> > Here what i done
> >
> >
> > close all;
> > clear all;
> > fo = 4; %frequency of the sine wave
> > Fs = 50; %sampling rate
> > Ts = 1/Fs; %sampling time interval
> > t = 0:Ts:1Ts; %sampling period
> > N=length(t);
> > y = 2*sin(2*pi*fo*t); %the sine curve
> >
> >
> > deltak=2*pi/N*deltax;
> > m=M/2:M/21
> > m*deltak=k (This is my frequency axis)
> > for i=1:N
> > X(i)=y(i)*exp(j*m*deltak*Ts)*Ts
> > end
> >
> >
> > I am doing some thing terribly wrong here ! I am confuse ...kindly help in writing code for DFT ...Assuming X as time and K as frequency in above expression......
> >
> > Thanks
> >
> > KK
>
> Hi KK, I see a couple problems with your DFT code. One problem is that in your for loop, you're not summing over the t variable. In the DFT at each frequency index, you have to sum over all values of the t variable, then you move to the next frequency index and sum over all values of t again. Below I calculate the DFT in a loop for your signal and then compare it against the builtin FFT routine in MATLAB.
>
> fo = 4; %frequency of the sine wave
> Fs = 50; %sampling rate
> Ts = 1/Fs; %sampling time interval
> t = 0:Ts:1Ts; %sampling period
> N=length(t);
> y = 2*sin(2*pi*fo*t); %the sine curve
> X = zeros(size(y));
> for k=1:N
> freq = (1j*2*pi*(k1))/N;
> Wnk = exp(freq.*(0:length(t)1));
> X(k)=Ts*sum(y.*Wnk);
> end
> Y = fft(y);
> subplot(211);
> plot(abs(X)); title('DFT');
> subplot(212);
> plot(abs(Y)); title('MATLAB FFT');
>
> Hope that helps,
> Wayne
