MATLAB Central Newsreader - ratio of fft's

ratio of fft's

Bruce — Tue, 20 May 2008 18:00:25 -0400

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

Re: ratio of fft's

roberson@ibd.nrc-cnrc.gc.ca (Walter Roberson) — Tue, 20 May 2008 18:54:28 -0400

In article <g0v3jp$27s$1@fred.mathworks.com>,
Bruce <italianasa84@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?

>% y2 = shifted sine wave
>y2 = sin(2*pi*10*t + phi).*hamming(length(t))';

When you use sin() to construct a sine wave that you submit to
fft, then due to round-off error in the calculation of the sine,
you will likely end up something that does not fft to a simple impulse
in the frequency domain: you will get a phase distortion that
will result in non-zero values at most points in the fft.
The power of that phase distortion is quite small, but it is present
and will have a noticable effect on your calculation.

If you start by -constructing- an impulse in the frequency domain
and ifft that back to a signal, the difference between the ifft'd
values and the sin() values you calculated will be on the order of
1E-14 to 1E-16. You cannot distinguish the two waveforms
if you plot them on top of each other: the changes are just too small.

The issue becomes quite prominent if you happen to use sind()
instead of sin(); some day, if I ever get in a fresh supply
of Round Tuits, I may investigate that further and possibly even
file a bug report; when I was experimenting with this a couple of
weeks ago, the sind() differences looked too large to be "reasonable".
--
"Tired minds don't plan well. Sleep first, plan later."
-- Walter Reisch

Re: ratio of fft's

dbd — Tue, 20 May 2008 21:38:32 -0400

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

Re: ratio of fft's

Ken Garrard — Wed, 21 May 2008 13:43:02 -0400

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

Re: ratio of fft's

Bruce — Thu, 22 May 2008 06:18:02 -0400

Thanks guys! That was exactly what I was looking for. If
anyone comes across this thread and wants more revised code
(there were just a couple small errors with the previous
post) my revised code is below. Thanks again!!

clear all; clc; close all;

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
nloop = 1000;
% 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 = cplxfilt(Y1);

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

    warning off MATLAB:divideByZero;
    ratio = unwrap(angle(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);
    plot(ratio*180/pi); %stem(ratio*180/pi);
    axis([-200, length(ratio)+200, -180, 180]);
    text(500,100,['degree shift:
',num2str(rem(phi*180/pi,360))]);
     text(500,50,['ratio: ',num2str(ratio(10)*180/pi)]);

    loops = loops + 1;
    pause(0.05);
end

function ff = cplxfilt(f)

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