Path: news.mathworks.com!not-for-mail
From: "Ken Garrard" <ken_garrardAT@ncsuDOT.edu>
Newsgroups: comp.soft-sys.matlab
Subject: Re: ratio of fft's
Date: Wed, 21 May 2008 13:43:02 +0000 (UTC)
Organization: North Carolina State University
Lines: 153
Message-ID: <g118t6$ll8$1@fred.mathworks.com>
References: <g0v3jp$27s$1@fred.mathworks.com> <b6c4d1df-b271-4336-a0e3-ba4e3f834342@l17g2000pri.googlegroups.com>
Reply-To: "Ken Garrard" <ken_garrardAT@ncsuDOT.edu>
NNTP-Posting-Host: webapp-03-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1211377382 22184 172.30.248.38 (21 May 2008 13:43:02 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Wed, 21 May 2008 13:43:02 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 216874
Xref: news.mathworks.com comp.soft-sys.matlab:469680


dbd <dbd@ieee.org> wrote in message <b6c4d1df-b271-4336-
a0e3-ba4e3f834342@l17g2000pri.googlegroups.com>...
> 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

Bruce,

In addition to properly using the window function as Dale 
suggests, you need to deal with the truncation effects of 
the FFT calculations before taking the ratio of the two 
complex vectors to construct a phase vector.  Otherwise the 
arc tangent of the ratio of two "noisy" values will produce 
the spikes that you observed.  Also do not force the FFT 
length to be a power of two.

Try the modified version of your code below and I think you 
will get what you expected; a spike at the appropriate 
frequency (positive and negative) and zeros everywhere else.
Without the hamming window the spike should be one sample 
wide.

Regards,
Ken



function phaseshift(nloop,usewin)

if nargin < 1 || isempty(nloop),   nloop = 1000; end
if nargin < 2 || isempty(usewin), usewin = true; end

Fs = 4000;      % sampling frequency
T = 1/Fs;       % sample time
L = 4000;       % length of signal
t = (0:L-1)*T;  % time vector

NFFT = L;       % length of FFT

% window function
if usewin, winy = hamming(length(t),'periodic').';
else       winy = ones(1,length(t));
end

% sin wave
y1 = sin(2*pi*10*t) .* winy;
Y1 = fft(y1,NFFT)/L;
Y1 = cplxf(Y1);

loops = 0;
while(loops < nloop);
    phi = loops*pi/360;
    y2 = sin(2*pi*10*t + phi) .* winy;
    Y2 = fft(y2,NFFT);
    Y2 = cplxf(Y2);

    warning off MATLAB:divideByZero;
    ratio = unwrap(angle(cplxfilt(Y2./Y1)));
    warning on MATLAB:divideByZero;

    subplot(1,2,1);
    plot(Fs*t,y1);     hold on;
    plot(Fs*t,y2,'g'); hold off;
    
    subplot(1,2,2);
    stem(ratio*180/pi);
    axis([-200, length(ratio)+200, -180, 180]);
    text(500,100,...
       ['degree shift: ',num2str(rem(phi*180/pi,360))]);
    
    loops = loops + 1;
    pause(.01);
end

function ff = cplxf(f,threshold)

error(nargchk(1,2,nargin));
if nargin < 2, threshold = 1e-10; end

% separate real and imaginary parts of input vector
fr = real(f);
fi = imag(f);

% find indices of values less than threshold
irz = find(abs(fr)<threshold);
iiz = find(abs(fi)<threshold);

% set "small" values to zero
fr(irz) = 0;
fi(iiz) = 0;

% rebuild complex vector
ff = fr + i*fi;
end

end