help-octave
[Top][All Lists]
Advanced

[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.



reply via email to

[Prev in Thread] Current Thread [Next in Thread]