From: kk KKsingh on
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,N-1} 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/2-1
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:1-Ts; %sampling period
N=length(t);
y = 2*sin(2*pi*fo*t); %the sine curve


deltak=2*pi/N*deltax;
m=-M/2:M/2-1
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
From: Rune Allnor on
On 20 Feb, 08:57, "kk KKsingh" <akikumar1...(a)gmail.com> wrote:
> I want to write my own code for DFT for regular sampling signall
....
> 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......

The error you make is that you think that a particular application
of the DFT requires a totally new implementation. It doesn't.

Use the standard (2D) DFT routines. Just add the book-keeping of the
frequency and wavenumber domains, which happen to be a couple of
vectors that have nothing to do with what goes on nside the (2D) DFT.

Rune
From: Wayne King on
"kk KKsingh" <akikumar1983(a)gmail.com> wrote in message <hlo4ke$4on$1(a)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,N-1} 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/2-1
> 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:1-Ts; %sampling period
> N=length(t);
> y = 2*sin(2*pi*fo*t); %the sine curve
>
>
> deltak=2*pi/N*deltax;
> m=-M/2:M/2-1
> 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 built-in FFT routine in MATLAB.

fo = 4; %frequency of the sine wave
Fs = 50; %sampling rate
Ts = 1/Fs; %sampling time interval
t = 0:Ts:1-Ts; %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*(k-1))/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
From: kk KKsingh on
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)=(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')




This is from Levent Sevgi
DOGUS University, Electronics and Communication Engineering Department
Zeamet Sokak, No 21, Ac&#305;badem, Istanbul, Turkey, lsevgi(a)dogus.edu.tr













"Wayne King" <wmkingty(a)gmail.com> wrote in message <hlolgh$eca$1(a)fred.mathworks.com>...
> "kk KKsingh" <akikumar1983(a)gmail.com> wrote in message <hlo4ke$4on$1(a)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,N-1} 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/2-1
> > 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:1-Ts; %sampling period
> > N=length(t);
> > y = 2*sin(2*pi*fo*t); %the sine curve
> >
> >
> > deltak=2*pi/N*deltax;
> > m=-M/2:M/2-1
> > 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 built-in FFT routine in MATLAB.
>
> fo = 4; %frequency of the sine wave
> Fs = 50; %sampling rate
> Ts = 1/Fs; %sampling time interval
> t = 0:Ts:1-Ts; %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*(k-1))/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
From: kk KKsingh on
Also I have few quick question assuming that my DFT expression looks like this


S(mdf)=sum(n=0 to N-1) s(ndt) exp(-j m df n dt) *dt


1. Above expression is for DFT
2. Now df=2*pi/N*dt

subtituting value from 2


S(mdf)=sum(n=0 to N-1) s(ndt) exp(j m (2 pi/(N dt)) n dt) *dt

So here is my Final expression for DFT

3. Now my question is that in the algorithm in most of the books i dont see this....infact what i see is very smiliar to you wrote...

freq = (-1j*2*pi*(k-1))/N;
Wnk = exp(freq.*(0:length(t)-1));
X(k)=Ts*sum(y.*Wnk);

why we are not using dt in denominater like in the above equation in this code. whats wrong in above equation. And when we take number of frequency points, i have seen some where that m=-M:M WHERE N=2M+1 so should i take K=1:M or K=-M:M .....i mean value of M comes in form of index as 1,2 3 or it comes like value at these index like m(1)=-50 m(2)=48. Actually I was trying to write code for irregular samples n now my mind is in dark black hole

"Wayne King" <wmkingty(a)gmail.com> wrote in message <hlolgh$eca$1(a)fred.mathworks.com>...
> "kk KKsingh" <akikumar1983(a)gmail.com> wrote in message <hlo4ke$4on$1(a)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,N-1} 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/2-1
> > 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:1-Ts; %sampling period
> > N=length(t);
> > y = 2*sin(2*pi*fo*t); %the sine curve
> >
> >
> > deltak=2*pi/N*deltax;
> > m=-M/2:M/2-1
> > 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 built-in FFT routine in MATLAB.
>
> fo = 4; %frequency of the sine wave
> Fs = 50; %sampling rate
> Ts = 1/Fs; %sampling time interval
> t = 0:Ts:1-Ts; %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*(k-1))/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