Prev: switching tv capture card pal<->ntsc standard
Next: Regarding xlswrite() function.
From: kk KKsingh on 1 Mar 2010 12:23
Since I am in FFT and DFT world now days..I was trying to read the code below few days back. According to the author x should be between -1/2 and 1/2..Just wondering if my t is between 0 to 2 then how can i get it between -1/2 and 1/2 without effecting the spectrum.....I hope i am clear in my question

And suppose if i run this code

clear all
close all
N=400; % Number sample points
fo = 10; %frequency of the sine wave
Fs = 100; %sampling rate
Ts = 1/Fs; %sampling time interval
t = 0:Ts:Ts*(N-1); %sampling period
y = 2*sin(2*pi*fo*t);

M=length(y);
N=M
freq=(-N/2):(N/2-1);

f_hat=y;

f=zeros(M,1);
f=exp(-2*pi*i*t'*freq)*f_hat';

stem(freqaxis1,abs(f))


I am not getting spectrum at right frequencies, I have taken this code from below

% NDFT.M
%
% Computes the nonequispaced dft and its adjoint.
%
% b=ndft(a,x,N,options,tflag)
%
% a Fourier coefficients or samples
% x nodes \subset [-1/2,1/2)
% N polynomial degree
% options .method 'direct' computes the columns of the nonequispaced
% Fourier matrix by direct calls of exp
% 'horner' avoids direct calls
% 'matrix' one matrix valued call of exp
%
% tflag 'notransp' evaluate trigonometric polynomial
% 'transp' adjoint transform
%
% Author: Stefan Kunis
%
% ---------------------------------------------------------------
%
% COPYRIGHT : (c) NUHAG, Dept.Math., University of Vienna, AUSTRIA
% http://nuhag.eu/
% Permission is granted to modify and re-distribute this
% code in any manner as long as this notice is preserved.
% All standard disclaimers apply.
%
function b=ndft(a,x,N,ft_opt,tflag)
%
% Computes the nonequispaced dft and its adjoint.
%
% b=ndft(a,x,N,options,tflag)
%
% a Fourier coefficients or samples
% x nodes \subset [-1/2,1/2)
% N polynomial degree
% options .method 'direct' computes the columns of the nonequispaced
% Fourier matrix by direct calls of exp
% 'horner' avoids direct calls
% 'matrix' one matrix valued call of exp
%
% tflag 'notransp' evaluate trigonometric polynomial
% 'transp' adjoint transform
%
% Author: Stefan Kunis

options=ft_opt.method;

M=length(x);
freq=(-N/2):(N/2-1);

if strcmp(tflag,'notransp')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% nonequispaced dft
%
% N/2-1
% f(j) = sum f_hat(k+N/2+1)*exp(-2*pi*i*k*x(j)), 1 <= j <= M.
% k=-N/2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

f_hat=a;
switch options
case 'direct'
f=zeros(M,1);
for k=1:N
f=f+f_hat(k).*exp(-2*pi*i*x*freq(k));
end;
case 'horner'
f=zeros(M,1);
exp_x=exp(2*pi*i*x);
for k=1:N
f=(f+f_hat(k)).*exp_x;
end;
f=f.*exp(-N*pi*i*x);
case 'matrix'
f=exp(-2*pi*i*x*freq)*f_hat;
otherwise
disp('unknown option')
end;
b=f;

elseif strcmp(tflag,'transp')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% adjoint nonequispaced dft
%
% M
% f_hat(k+N/2+1) = sum f(j)*exp(2*pi*i*k*x(j)), -N/2 <= k <= N/2-1.
% j=1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

f=a;
switch options
case 'direct'
f_hat=zeros(N,1);
for k=1:N
f_hat(k)=exp(-2*pi*i*x*freq(k))' *f;
% exp(-2*pi*i*freq(k)*x)' *f differs by 10^-11!?
end;
case 'horner'
f_hat=zeros(N,1);
f=f.*exp(-N*pi*i*x);
exp_x=exp(2*pi*i*x);
for k=1:N
f_hat(k)=sum(f);
f=f.*exp_x;
end;
case 'matrix'
f_hat=(exp(-2*pi*i*x*freq)') *f;
otherwise
disp('unknown option')
end;
b=f_hat;

end;
From: Greg Heath on 4 Mar 2010 03:31
On Mar 1, 12:23 pm, "kk KKsingh" <akikumar1...(a)gmail.com> wrote:
-----SNIP
clear all
close all
clc
N = 40; % Number sample points
fo = 10; %frequency of the sine wave
Fs = 100; %sampling rate
Ts = 1/Fs; %sampling time interval
t = 0:Ts:Ts*(N-1); %sampling period
y = 2*sin(2*pi*fo*t);

figure(1)
plot(t,y,'--.')

freq1 = (Fs/N)*( 0 : N-1 );
freq2 = (Fs/N)*(-N/2 : N/2-1);

Y1 = exp(-2*pi*i*t'*freq1)*y'/N;
Y2 = fftshift(Y1);

Y3 = exp(-2*pi*i*t'*freq2)*y'/N;

' NOTE 1 : Y3 ~= Y2 (See Below) '

figure(2)
subplot(2,1,1)
stem(freq1,abs(Y1))
subplot(2,1,2)
stem(freq2,abs(Y2))

' NOTE 2 '

maxabs(Y3(1:2:end) - Y1(1:2:end)) % 1.8262e-015
maxabs(Y3(2:2:end) + Y1(2:2:end)) % 1.6941e-015

' NOTE 3 '
' What happens to the spectrum when N < 40 ?) '


Hope this helps.

Greg
 | 
Pages: 1
Prev: switching tv capture card pal<->ntsc standard
Next: Regarding xlswrite() function.