Path: news.mathworks.com!newsfeed-00.mathworks.com!newsfeed2.dallas1.level3.net!news.level3.com!postnews.google.com!l17g2000pri.googlegroups.com!not-for-mail
From: dbd <dbd@ieee.org>
Newsgroups: comp.soft-sys.matlab
Subject: Re: ratio of fft's
Date: Tue, 20 May 2008 14:38:32 -0700 (PDT)
Organization: http://groups.google.com
Lines: 52
Message-ID: <b6c4d1df-b271-4336-a0e3-ba4e3f834342@l17g2000pri.googlegroups.com>
References: <g0v3jp$27s$1@fred.mathworks.com>
NNTP-Posting-Host: 75.17.140.180
Mime-Version: 1.0
Content-Type: text/plain; charset=ISO-8859-1
Content-Transfer-Encoding: 7bit
X-Trace: posting.google.com 1211319513 19157 127.0.0.1 (20 May 2008 21:38:33 GMT)
X-Complaints-To: groups-abuse@google.com
NNTP-Posting-Date: Tue, 20 May 2008 21:38:33 +0000 (UTC)
Complaints-To: groups-abuse@google.com
Injection-Info: l17g2000pri.googlegroups.com; posting-host=75.17.140.180; 
User-Agent: G2/1.0
X-HTTP-UserAgent: Mozilla/5.0 (Windows; U; Windows NT 5.1; en-US; rv:1.8.1.14) 
Xref: news.mathworks.com comp.soft-sys.matlab:469558


On May 20, 11:00 am, "Bruce " <italianas...@gmail.com> wrote:
> A simple divison, on paper, of two fft's proves that the
> ratio of the fourier transform of two sine waves with the
> same frequency gives exp(-j*phi) where phi is the phase
> shift between the two. Therefore if you take
> angle(exp(-j*phi)) you should get phi, which is a constant.
> Why in matlab do you not get a constant?? Where do those
> spikes come from?
>
> % start of code. press spacebar to continue shifting.
> clear all; close all; clc;
> Fs = 4000; % sampling frequency
> T = 1/Fs;  % sample time
> L = 4000; %length of signal
> t = (0:L-1)*T;  % time vector
> NFFT = 2^nextpow2(L); % Next power of 2 from length of y
>
> y1 =  sin(2*pi*10*t).* hamming(length(t))';
> Y1 = fft(y1,NFFT)/L;
> f = Fs/2*linspace(0,1,NFFT/2);
> loops = 1;
> phi = 0;
> % y2 = shifted sine wave
> y2 =  sin(2*pi*10*t + phi).*hamming(length(t))';
> Y2 = fft(y2,NFFT);
>
> while(loops < 10000);
>     phi = loops*pi/360;
>     y2 =  sin(2*pi*10*t + phi) .* hamming(length(t))';
>     Y2 = fft(y2,NFFT);
>
>     subplot(1,2,1);
>     plot(Fs*t,y1); hold on; plot(Fs*t,y2,'g'); hold off;
>
>     subplot(1,2,2);    ratio = unwrap(angle(Y2./Y1));
>     plot(ratio); axis([0, length(ratio), -5, 5]);
>     text(500,4, ['degree shift: ',
> num2str(rem(rad2deg(phi),360))]);
>
>     loops = loops + 1;
>     pause();
> end

If you are going to window the data you must window it properly to get
proper results. That is to get bin centered tones to be bin centered,
etc.

In the hamming calls use the SFLAG parameter set to 'periodic'. See
the hamming help. The default is not correct for windowing.

Dale B. Dalrymple
http://dbdimages.com