[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Problem with using fft and ifft functions
From: |
babelproofreader |
Subject: |
Problem with using fft and ifft functions |
Date: |
Tue, 21 Jul 2009 07:37:15 -0700 (PDT) |
I'm trying to recreate code in Octave from a C++ file to implement a lowpass
smoothing filter which uses fast fourier transforms and inverse fast fourier
transforms. Now admittedly my knowledge of DSP is very basic, but I
understood that a fast fourier transform transforms a signal from the time
domain into the frequency domain and the inverse fast fourier transform
transforms it back from the frequency domain to the time domain. However,
when I use the functions fft and then ifft something seems to go wrong, but
in my ignorance I can't figure out how to fix it, or even what is going
wrong. I also have to say that my knowledge of C++ is almost non existant.
The section of relevant C++ code is
double slope = 0; // will be modified on call to detrend
double intercept = 0;
int length = 0; // original caller size
int n = 0; // size raised to next power of 2 for fft
int i = 0;
length = in->getSize();
// Detrend input series
PlotLine *series = detrend(in, slope, intercept, true);
// Raise length to next power of 2, pad with zero
PlotLine *series2 = raise2Power(series, 0);
n = series2->getSize();
//qtsFFT fft(n); // construct fft object
fft = new qtsFFT(n);
// do fft
PlotLine * fftFreq = fft->do_FFTqts(series2);
//PlotLine * fftFreq = fft.do_FFTqts(series2);
// apply low pass filter
double f = 0;
double dist = 0;
double wt = 0;
int halfn = n/2;
double freqSave = fftFreq->getData(halfn);
for (i = 0 ; i < halfn ; i++)
{
f = (double) i / (double) n ; // Frequency
if (f <= fre) // Flat response
wt = 1.0 ;
else
{
dist = (f - fre) / wid;
wt = exp ( -dist * dist ) ;
}
fftFreq->setData(i, fftFreq->getData(i) * wt) ;
fftFreq->setData(halfn + i, fftFreq->getData(halfn + i) * wt) ;
}
dist = (0.5 - fre) / wid; // Do Nyquist in fftFreq[0]
fftFreq->setData(halfn, freqSave * exp ( -dist * dist )) ;
// Do inverse FFT to recover real domain
PlotLine *fftReal = fft->do_iFFTqts(fftFreq);
//PlotLine *fftReal = fft.do_iFFTqts(fftFreq);
// Retrend input series, n.b. original length
PlotLine *series3 = detrend(fftReal, slope, intercept, false);
for (i = 0; i < length; i++)
out->append(series3->getData(i));
delete series;
delete series2;
delete series3;
delete fftReal;
delete fftFreq;
delete fft;
return out;
}
and my messy working file attempt at an Octave implementation, stopping at
the ifft failure is
data=load("-ascii","gc");
n=rows(data);
dn=datenum(data(:,1),data(:,2),data(:,3)); % col 1 is year, col 2 is month,
col 1 is day
dy=data(:,3);
mth=data(:,2);
yr=data(:,1);
open=data(:,4);
high=data(:,5);
low=data(:,6);
close=data(:,7);
typprice=(high.+low.+close)./3;
fre=0.05;
wid=0.2;
x=(1:n)';
p = polyfit(x,typprice,1);
lmsfit=polyval(p,x);
lmstypprice=typprice.-lmsfit;
if rem(n,2)==0
fn=n;
else
fn=n+1;
endif
[a,b]=stft(lmstypprice);
fftFreq=zeros(fn,1);
%plot(b);
halfn=fn/2;
freqsave=a(halfn,1);
% for i=1:halfn-1
% f=i/fn; % frequency
% if (f<fre)==1 % flat response set at 0.05
% wt=1;
% else
% dist=(f-fre)/wid;
% wt=exp(-dist*dist);
% endif
% fftFreq(i,1)=fftFreq(i,1)*wt;
% fftFreq(i+halfn,1)=fftFreq(i+halfn,1)*wt;
% endfor
for i=1:halfn-1
f=i/fn; % frequency
if (f<fre)==1 % flat response set at 0.05
wt=1;
else
dist=(f-fre)/wid;
wt=exp(-dist*dist);
endif
fftFreq(i,1)=a(i,1)*wt;
fftFreq(i+halfn,1)=a(i+halfn,1)*wt;
endfor
dist=(0.5-fre)/wid; % Do Nyquist in fftFreq[0]
fftFreq(halfn,1)=freqsave*exp(-dist*dist);
remade=zeros(fn,1);
%plot(fftFreq);
plot(remade);
% Do inverse FFT to recover real domain
fftReal=ifft(fftFreq);
plot(fftReal);
A=[x,open,high,low,close,lmsfit];
dlmwrite("chart",A)
If anyone could give me advice, I'd really appreciate it.
--
View this message in context:
http://www.nabble.com/Problem-with-using-fft-and-ifft-functions-tp24588890p24588890.html
Sent from the Octave - General mailing list archive at Nabble.com.
- Problem with using fft and ifft functions,
babelproofreader <=