Path: news.mathworks.com!not-for-mail
From: "Bruce " <italianasa84@gmail.com>
Newsgroups: comp.soft-sys.matlab
Subject: Re: ratio of fft's
Date: Thu, 22 May 2008 06:18:02 +0000 (UTC)
Organization: NASA
Lines: 73
Message-ID: <g1336q$31r$1@fred.mathworks.com>
References: <g0v3jp$27s$1@fred.mathworks.com> <b6c4d1df-b271-4336-a0e3-ba4e3f834342@l17g2000pri.googlegroups.com> <g118t6$ll8$1@fred.mathworks.com>
Reply-To: "Bruce " <italianasa84@gmail.com>
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 1211437082 3131 172.30.248.38 (22 May 2008 06:18:02 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Thu, 22 May 2008 06:18:02 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 774179
Xref: news.mathworks.com comp.soft-sys.matlab:469800


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