NDFT
in [Matlab]
|
Prev: Bisection Method
Next: GUI with wavelet transformates
From: Greg Heath on 1 Jun 2010 00:28 >Newsgroups: comp.soft-sys.matlab >From: "kk KKsingh" <akikumar1...(a)gmail.com> >Date: Mon, 1 Mar 2010 17:23:22 +0000 (UTC) >Local: Mon, Mar 1 2010 1:23 pm >Subject: NDFT > >Since I am in FFT and DFT world now days... >-----SNIP >And suppose if i run this code >-----SNIP >I am not getting spectrum at right frequencies,... Try this: close all, clear all, clc, k=0 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); k=k+1, figure(k) % Figure 1 plot(t,y) title('Original Time Signal') % Standard Unipolar DFT Freq Range freq1 = (Fs/N)*( 0 : N-1 ); % Standard Bipolar DFT Freq Range freq2 = (Fs/N)*(-N/2 : N/2-1); % Y1 defined on freq1 Y1 = exp(-2*pi*i*t'*freq1)*y'/N; % Y2 defined on freq2 Y2 = fftshift(Y1); % Y3 also assumed to be defined on freq2 Y3 = exp(-2*pi*i*t'*freq2)*y'/N; ' HOWEVER! ... ' A1 = abs(Y1) ; A2 = abs(Y2) ; A3 = abs(Y3) ; k=k+1, figure(k) % Figure 2 subplot(3,1,1) stem(freq1,A1) ylabel('abs(Y1)') title('Magnitude vs Frequency') subplot(3,1,2) stem(freq2,A2,'g') ylabel('abs(Y2)') subplot(3,1,3) stem(freq2,A3,'r') ylabel('abs(Y3)') disp([... 'Therefore, ';... ' 1. Y3 ~= Y2 ';... ' 2. It looks like ';... ' a. Y3 should be defined on freq1! ';... ' b. abs(Y3) = abs(Y1) ';... ' c. angle(Y3) and angle(Y1) ';... ' should have a deterministic ';... ' relationship ';... ' ';... ' Checking: ']) 'Find Peak Locations' [A1max ipeak1] = max(A1) % [ 1 37 ] ipeaks1 = find(A1max-A1 < 1e-15) % [ 5 37]' inonpeak1 = setdiff(1:N,ipeaks1); [A3max ipeak3] = max(A3) % [ 1 37 ] ipeaks3 = find(A3max-A3 < 1e-15) % [ 5 37]' inonpeak3 = setdiff(1:N,ipeaks3); R1 = real(Y1); I1 = imag(Y1); P1 = angle(Y1); R3 = real(Y3); I3 = imag(Y3); P3 = angle(Y3); checkY = maxabs(Y3-Y1) % 5.459e-015 checkR = maxabs(R3-R1) % 1.8074e-015 checkI = maxabs(I3-I1) % 5.3318e-015 checkA = maxabs(A3-A1) % 1.6335e-015 checkP = maxabs(P3-P1) % 4.8608 'Still can''t discern the relationship between P3 and P1!' k=k+1, figure(k) % Figure 3 subplot(221), hold on plot(freq1,R1) plot(freq1,R3,'r') legend('Y1','Y3',4) title('Real Parts') subplot(222), hold on plot(freq1,I1) plot(freq1,I3,'r') title('Imaginary Parts') subplot(223), hold on plot(freq1,A1) plot(freq1,A3,'r') title('Magnitudes') subplot(224), hold on plot(freq1,P1) plot(freq1,P3,'r') title('Phases') ' AHA! Note that' maxabsR1 = maxabs(R1) % 2.0264e-015 maxabsR3 = maxabs(R3) % 3.6152e-016 maxabsI1nonpeaks = maxabs(I1(inonpeak1)) % 2.5174e-015 maxabsI3nonpeaks = maxabs(I3(inonpeak3)) % 2.8144e-015 disp('Therefore the nonpeak phase estimates are unreliable!') disp('Correcting ... ') R1 = zeros(N,1); I1(inonpeak1) = zeros(N-2,1); R3 = zeros(N,1); I3(inonpeak1) = zeros(N-2,1); Y1 = complex(R1,I1); Y3 = complex(R3,I3); P1 = angle(Y1); P3 = angle(Y3); checkY = maxabs(Y3-Y1) % 5.3318e-015 checkR = maxabs(R3-R1) % 0 checkI = maxabs(I3-I1) % 5.5511e-015 checkA = maxabs(A3-A1) % 1.6335e-015 checkP = maxabs(P3-P1) % 0 k=k+1, figure(k) % Figure 4 subplot(221), hold on plot(freq1,R1) plot(freq1,R3,'r') legend('Y1','Y3',4) title('Corrected Real Parts') subplot(222), hold on plot(freq1,I1) plot(freq1,I3,'r') title('Corrected Imaginary Parts') subplot(223), hold on plot(freq1,A1) plot(freq1,A3,'r') title('Magnitudes') subplot(224), hold on plot(freq1,P1) plot(freq1,P3,'r') title('Phases') U1 = unwrap(P1); U3 = unwrap(P3); U31 = unwrap(P3-P1); k=k+1, figure(k) % Figure 5 subplot(221), hold on plot(freq1,P1) plot(freq1,P3,'r') title('Phase vs Frequency') subplot(222), hold on plot(freq1,U1) plot(freq1,U3,'r') title('Unwrapped Phase vs Frequency') subplot(223), hold on plot(freq1,P3-P1) title('Phase Difference vs Frequency') subplot(224), hold on plot(freq1,U31) title('Unwrapped Phase Difference vs Frequency') sortfigs(1,3) Hope this helps. Greg
From: Greg Heath on 1 Jun 2010 00:33 On Jun 1, 12:28 am, Greg Heath <he...(a)alumni.brown.edu> wrote: > sortfigs(1,3) Correction: sortfigs(1,5) function sortfigs(kmin,kmax) for j = kmax:-1:kmin figure(j) end Hope this helps. Greg
From: Greg Heath on 1 Jun 2010 10:50 On Jun 1, 12:28 am, Greg Heath <he...(a)alumni.brown.edu> wrote: > >Newsgroups: comp.soft-sys.matlab > >From: "kk KKsingh" <akikumar1...(a)gmail.com> > >Date: Mon, 1 Mar 2010 17:23:22 +0000 (UTC) > >Local: Mon, Mar 1 2010 1:23 pm > >Subject: NDFT I forgot to mention that because of the length of time since the original NDFT March 1 thread http://groups.google.com/group/comp.soft-sys.matlab/ browse_thread/thread/8af150f9618bd587/d509f46d08ccee7e I had to start this new thread with the same name > >Since I am in FFT and DFT world now days... > >-----SNIP > >And suppose if i run this code > >-----SNIP > >I am not getting spectrum at right frequencies,... > > Try this: > close all, clear all, clc, k=0 > > 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); > > k=k+1, figure(k) % Figure 1 > plot(t,y) > title('Original Time Signal') > > % Standard Unipolar DFT Freq Range > freq1 = (Fs/N)*( 0 : N-1 ); > % Standard Bipolar DFT Freq Range > freq2 = (Fs/N)*(-N/2 : N/2-1); > > % Y1 defined on freq1 > Y1 = exp(-2*pi*i*t'*freq1)*y'/N; > % Y2 defined on freq2 > Y2 = fftshift(Y1); > % Y3 also assumed to be defined on freq2 > Y3 = exp(-2*pi*i*t'*freq2)*y'/N; > > ' HOWEVER! ... ' -----SNIP > ' 1. Y3 ~= Y2 ';... > ' 2. It looks like ';... > ' a. Y3 should be defined on freq1! ';... > ' b. abs(Y3) = abs(Y1) ';... > ' c. angle(Y3) and angle(Y1) ';... > ' should have a deterministic ';... > ' relationship ';... > ' ';... > ' Checking: ']) -----SNIP Bottom line: Y3 = Y1 Question: WHY??? Answer: Because the OP used row vectors instead of column vectors which, apparently led him to use the incorrect formulation W = exp( -2*pi*i * t' * freq ) Instead of the correct formulation W = exp( -2*pi*i * freq' * t) Making the correction yields Y3 = Y2. Hope this helps. Greg
|
Pages: 1 Prev: Bisection Method Next: GUI with wavelet transformates |