I'm still having problems. My .oct function code now is
#include <octave/oct.h>
#include <octave/dColVector.h>
#include <octave/parse.h>
#include <math.h>
DEFUN_DLD (pass4, args, , "Help String")
{
octave_value retval;
ColumnVector detrend_input = args(0).column_vector_value ();
ColumnVector detrend_output(detrend_input);
double intercept = detrend_input(0);
double slope = ( detrend_input( detrend_input.length () - 1 ) - intercept
) / ( detrend_input.length () - 1 );
octave_value_list fft_result;
ComplexColumnVector freq_rep;
octave_value_list ifft_result;
octave_value_list real_result;
ColumnVector real_result_rep;
ColumnVector retrended_real_result(detrend_input);
for (octave_idx_type ii (0); ii<detrend_input.length (); ii++)
{
detrend_output(ii) = detrend_input(ii) - intercept - (ii) * slope;
}
fft_result = feval ("fft", octave_value (detrend_output), 1);
freq_rep = fft_result (0).complex_column_vector_value ();
// The filter code
ComplexColumnVector filtered_freq_rep(freq_rep);
double freq = 0.05;
double wdth = 0.2;
double f = 0.0;
double wt = 0.0;
double dist = 0.0;
for (octave_idx_type ii (1); ii < (detrend_input.length () / 2 ); ii++)
{
f = (ii) / detrend_input.length ();
if (f <= freq)
{
wt = 1;
}
else
{
dist = ( f - freq ) / wdth;
wt = exp ( -dist * dist );
}
filtered_freq_rep(ii) = freq_rep(ii) * wt;
filtered_freq_rep(detrend_input.length () - ii ) =
freq_rep(detrend_input.length () - ii ) * wt;
}
dist = ( 0.5 - freq ) / wdth; // Do Nyquist
filtered_freq_rep(detrend_input.length () / 2 ) =
freq_rep(detrend_input.length () / 2 ) * exp ( -dist * dist);
// end of filer code
ifft_result = feval ("ifft", octave_value (filtered_freq_rep), 1);
real_result = feval ("real", ifft_result, 1);
real_result_rep = real_result (0).column_vector_value ();
for (octave_idx_type ii (0); ii<detrend_input.length (); ii++)
{
retrended_real_result(ii) = real_result_rep(ii) + intercept + (ii) *
slope;
}
retval = retrended_real_result;
return retval;
}
which compiles without giving any error messages. However, when I run it
with data the output is exactly the same as the input - no filtering appears
to be taking place. The code in m that I'm trying to implement is
function y=fftlowpass(x,fre,wid)
%function y=fftlowpass(x,fre,wid)
if(exist('wid')~=1) wid=0.2; end
if(exist('fre')~=1) fre=0.05; end
L=length(x);
[y,slope,intercept]=detrend(x,0,0,1);
L2=2^nextpow2(L); half_L2=L2/2;
y2=zeros(1,L2); y2(1:L)=y;
Y2=fft(y2);
for(i=2:half_L2)
f=(i-1)/L2;
if(f<=fre)
wt=1;
else
dist=(f-fre)/wid;
wt=exp(-dist*dist);
end
Y2(i)=Y2(i)*wt;
Y2(L2-i+2)=Y2(L2-i+2)*wt;
end
dist=(0.5-fre)/wid;
Y2(half_L2+1)=Y2(half_L2+1)*exp(-dist*dist);
y2=real(ifft(Y2));
y=detrend(y2(1:L),slope,intercept,0);
function [y,slope,intercept]=detrend(x,slope,intercept,flag)
L=length(x);
y=zeros(1,L);
if(flag==1)
intercept=x(1);
slope=(x(L)-intercept)/(L-1);
for(i=1:L)
y(i)=x(i)-intercept-(i-1)*slope;
end
else
for(i=1:L)
y(i)=x(i)+intercept+(i-1)*slope;
end
end
I have not bothered to implement that part of the code which raises to the
next power of 2 as I know that my data will not need this. I have checked
and as far as I can see my .oct code is a faithful translation of the m
code, so I am perplexed as to why the output results are not the same for
each.