How do you pass two Octave functions to an .oct function?

[babelproofreader]

I would like to pass two Octave functions, specifically fft and ifft, to an .oct file as part of a compiled function I wish to write. I have read the relevant on-line manual appendix and know the syntax to pass one Octave function to an .oct file, but how do I pass two functions?

[Jaroslav Hajek-2]

On Tue, Oct 13, 2009 at 12:53 PM, babelproofreader
<babelproofreader@...> wrote:
>
> I would like to pass two Octave functions, specifically fft and ifft, to an
> .oct file as part of a compiled function I wish to write. I have read the
> relevant on-line manual appendix and know the syntax to pass one Octave
> function to an .oct file, but how do I pass two functions?

Please be more specific and show some code if possible. You can of
course pass any number of function handles to a compiled functions,
either as individual arguments or, for example, in a cell array.

regards

--
RNDr. Jaroslav Hajek
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[babelproofreader]

>Please be more specific and show some code if possible...

My code so far is

#include <octave/oct.h>
#include <octave/dColVector.h>
#include <octave/parse.h>
     
DEFUN_DLD (rollingfft, args, , "Help String")
{
octave_value_list retval;
octave_value_list fft_result_vector;
octave_value_list ifft_result_vector;
octave_value_list real_result_vector;

 octave_function *fcn1 = args(0).function_value (); // supply Octave "fft" function
 octave_function *fcn2 = args(1).function_value (); // supply Octave "ifft" function
 octave_function *fcn3 = args(2).function_value (); // supply Octave "real" function
 ColumnVector input_vector = args(3).column_vector_value ();
 ColumnVector output_vector = args(3).column_vector_value ();
 OCTAVE_LOCAL_BUFFER( double , rolling_vector , 32 );

   for (octave_idx_type ii (31); ii<input_vector.length (); ii++)                                        
     {
       for (octave_idx_type jj (0); jj<(32); jj++)
         {
         rolling_vector(jj) = input_vector(ii - (ii - jj));
         }
           // detrend rolling vector code here

           fft_result_vector = feval( fcn1 , rolling_vector );
 
           // fft_result_vector frequency filter code here

           ifft_result_vector = feval( fcn2 , fft_result_vector );
           real_result_vector = feval( fcn3 , ifft_result_vector );

           // real_result_vector retrend code here
       
           output_vector(ii) = real_result_vector(31);
     }
                                   

 retval = output_vector;                                                                        

return retval;                                                                      
}

but it fails on compile with

octave:1> mkoctfile rollingfft.cc
rollingfft.cc: In function ‘octave_value_list Frollingfft(const octave_value_list&, int)’:
rollingfft.cc:27: error: call of overloaded ‘feval(octave_function*&, double*&)’ is ambiguous
/usr/include/octave-3.0.1/octave/parse.h:123: note: candidates are: octave_value_list feval(octave_function*, const octave_value_list&, int) <near match>
/usr/include/octave-3.0.1/octave/parse.h:126: note:                 octave_value_list feval(const octave_value_list&, int) <near match>
rollingfft.cc:36: error: cannot convert ‘octave_value’ to ‘double’ in assignment
rollingfft.cc:40: error: no match for ‘operator=’ in ‘retval = output_vector’
/usr/include/octave-3.0.1/octave/oct-obj.h:75: note: candidates are: octave_value_list& octave_value_list::operator=(const octave_value_list&)

What I am attempting to do is to smooth a time series by moving a window across the series, and using fft to lowpass filter in the frequency domain.

[Jaroslav Hajek-2]

On Thu, Oct 15, 2009 at 3:40 PM, babelproofreader
<babelproofreader@...> wrote:

>
>>Please be more specific and show some code if possible...
>
> My code so far is
>
> #include <octave/oct.h>
> #include <octave/dColVector.h>
> #include <octave/parse.h>
>
> DEFUN_DLD (rollingfft, args, , "Help String")
> {
> octave_value_list retval;
> octave_value_list fft_result_vector;
> octave_value_list ifft_result_vector;
> octave_value_list real_result_vector;
>
>  octave_function *fcn1 = args(0).function_value (); // supply Octave "fft"
> function
>  octave_function *fcn2 = args(1).function_value (); // supply Octave "ifft"
> function
>  octave_function *fcn3 = args(2).function_value (); // supply Octave "real"
> function
>  ColumnVector input_vector = args(3).column_vector_value ();
>  ColumnVector output_vector = args(3).column_vector_value ();
>  OCTAVE_LOCAL_BUFFER( double , rolling_vector , 32 );
>
>   for (octave_idx_type ii (31); ii<input_vector.length (); ii++)
>     {
>       for (octave_idx_type jj (0); jj<(32); jj++)
>         {
>         rolling_vector(jj) = input_vector(ii - (ii - jj));
>         }
>           // detrend rolling vector code here
>
>           fft_result_vector = feval( fcn1 , rolling_vector );
>
>           // fft_result_vector frequency filter code here
>
>           ifft_result_vector = feval( fcn2 , fft_result_vector );
>           real_result_vector = feval( fcn3 , ifft_result_vector );
>
>           // real_result_vector retrend code here
>
>           output_vector(ii) = real_result_vector(31);
>     }
>
>
>  retval = output_vector;
>
> return retval;
> }
>
> but it fails on compile with
>
> octave:1> mkoctfile rollingfft.cc
> rollingfft.cc: In function ‘octave_value_list Frollingfft(const
> octave_value_list&, int)’:
> rollingfft.cc:27: error: call of overloaded ‘feval(octave_function*&,
> double*&)’ is ambiguous
> /usr/include/octave-3.0.1/octave/parse.h:123: note: candidates are:
> octave_value_list feval(octave_function*, const octave_value_list&, int)
> <near match>
> /usr/include/octave-3.0.1/octave/parse.h:126: note:
> octave_value_list feval(const octave_value_list&, int) <near match>
> rollingfft.cc:36: error: cannot convert ‘octave_value’ to ‘double’ in
> assignment
> rollingfft.cc:40: error: no match for ‘operator=’ in ‘retval =
> output_vector’
> /usr/include/octave-3.0.1/octave/oct-obj.h:75: note: candidates are:
> octave_value_list& octave_value_list::operator=(const octave_value_list&)
>

I'd think these messages are pretty understandable.
You simply can't pass a raw double * pointer (the local buffer) to
feval, because it can't handle it. You need to create a valid
octave_value from it. Declare rolling_vector as ColumnVector and it
will likely work.

> What I am attempting to do is to smooth a time series by moving a window
> across the series, and using fft to lowpass filter in the frequency domain.

Have you tried to do this in m-code first? Maybe it could be made fast
enough. Especially if you upgrade Octave (3.0.1 is ancient).

hth

--
RNDr. Jaroslav Hajek
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[babelproofreader]

>I'd think these messages are pretty understandable.
>You simply can't pass a raw double * pointer (the local buffer) to
>feval, because it can't handle it. You need to create a valid
>octave_value from it. Declare rolling_vector as ColumnVector and it
>will likely work.

I have stripped my code down to the bare minimum, declared the input as ColumnVector, and it still fails on compile. My revised code is

#include <octave/oct.h>
#include <octave/dColVector.h>
#include <octave/parse.h>
     
DEFUN_DLD (rollingfft, args, , "Help String")
{
octave_value_list retval;
octave_value_list fft_result_vector;
octave_value_list ifft_result_vector;
octave_value_list real_result_vector;

 octave_function *fcn1 = args(0).function_value ();                                // supply Octave "fft" function
 octave_function *fcn2 = args(1).function_value ();                                // supply Octave "ifft" function
 octave_function *fcn3 = args(2).function_value ();                                // supply Octave "real" function
 ColumnVector input_vector = args(3).column_vector_value ();
 ColumnVector output_vector = args(3).column_vector_value ();

 fft_result_vector = feval( fcn1 , input_vector );
 ifft_result_vector = feval( fcn2 , fft_result_vector );
 real_result_vector = feval( fcn3 , ifft_result_vector );
 retval = real_result_vector;                                                                        

return retval;                                                                      
}

and the error message is

octave:1> mkoctfile rollingfft.cc
rollingfft.cc: In function ‘octave_value_list Frollingfft(const octave_value_list&, int)’:
rollingfft.cc:18: error: no matching function for call to ‘feval(octave_function*&, ColumnVector&)’
/usr/include/octave-3.0.1/octave/parse.h:118: note: candidates are: octave_value_list feval(const std::string&, const octave_value_list&, int)
/usr/include/octave-3.0.1/octave/parse.h:123: note:                 octave_value_list feval(octave_function*, const octave_value_list&, int)
/usr/include/octave-3.0.1/octave/parse.h:126: note:                 octave_value_list feval(const octave_value_list&, int)
 

[David Grundberg-2]

babelproofreader skrev:

>> I'd think these messages are pretty understandable.
>> You simply can't pass a raw double * pointer (the local buffer) to
>> feval, because it can't handle it. You need to create a valid
>> octave_value from it. Declare rolling_vector as ColumnVector and it
>> will likely work.
>>    
>
> I have stripped my code down to the bare minimum, declared the input as
> ColumnVector, and it still fails on compile. My revised code is
>
> #include <octave/oct.h>
> #include <octave/dColVector.h>
> #include <octave/parse.h>
>      
> DEFUN_DLD (rollingfft, args, , "Help String")
> {
> octave_value_list retval;
> octave_value_list fft_result_vector;
> octave_value_list ifft_result_vector;
> octave_value_list real_result_vector;
>
>  octave_function *fcn1 = args(0).function_value ();                              
> // supply Octave "fft" function
>  octave_function *fcn2 = args(1).function_value ();                              
> // supply Octave "ifft" function
>  octave_function *fcn3 = args(2).function_value ();                              
> // supply Octave "real" function
>  ColumnVector input_vector = args(3).column_vector_value ();
>  ColumnVector output_vector = args(3).column_vector_value ();
>
>  fft_result_vector = feval( fcn1 , input_vector );
>  ifft_result_vector = feval( fcn2 , fft_result_vector );
>  real_result_vector = feval( fcn3 , ifft_result_vector );
>  retval = real_result_vector;                                                                        
>
> return retval;                                                                      
> }
>
> and the error message is
>
> octave:1> mkoctfile rollingfft.cc
> rollingfft.cc: In function ‘octave_value_list Frollingfft(const
> octave_value_list&, int)’:
> rollingfft.cc:18: error: no matching function for call to
> ‘feval(octave_function*&, ColumnVector&)’
> /usr/include/octave-3.0.1/octave/parse.h:118: note: candidates are:
> octave_value_list feval(const std::string&, const octave_value_list&, int)
> /usr/include/octave-3.0.1/octave/parse.h:123: note:                
> octave_value_list feval(octave_function*, const octave_value_list&, int)
> /usr/include/octave-3.0.1/octave/parse.h:126: note:                
> octave_value_list feval(const octave_value_list&, int)
>  
>  

I don't understand under which circumstances you'd want to replace the
fourier transform functions. Passing the functions like arguments makes
me suspect you haven't looked at the API. Also, naming a variable of
type octave_value_list "vector" is confusing, so I removed that.

Have you tried implementing this in Octave code (.m)?


#include <octave/oct.h>
#include <octave/dColVector.h>
#include <octave/parse.h>

DEFUN_DLD (rollingfft, args, , "Help String")
{
  octave_value_list retval;
  octave_value_list fft_result;
  ComplexColumnVector freq_rep;
  octave_value_list ifft_result;
  octave_value_list real_result;

  fft_result = feval ("fft", args (0), 1);
  freq_rep = fft_result (0).complex_column_vector_value ();
  ifft_result = feval ("ifft", octave_value (freq_rep), 1);
  real_result = feval ("real", ifft_result, 1);
  retval = real_result;

  return retval;
}



_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[babelproofreader]

>Have you tried implementing this in Octave code (.m)?

Yes, I have; the code is below.

function smoothed_output=fastrollingfftlowpass(data,fre,wid)
%function y=fftlowpass(rolling_vector,fre,wid)
if(exist('wid')~=1) wid=0.2; endif
if(exist('fre')~=1) fre=0.05; endif
n=rows(data);
smoothed_output=data;
L=16;
L_2=L/2;
rolling_vector=zeros(L,1);
detrend_rolling_vector=zeros(L,1);

for (ii=L:n)

rolling_vector(1:L,1)=data(ii-(L-1):ii,1);

% detrend the rolling_vector %
intercept=rolling_vector(1,1);
slope=(rolling_vector(L,1)-intercept)/(L-1);
detrend_rolling_vector=detrend_vector(rolling_vector); % calls an .oct file "detrend_vector"
% end of detrend code %

% the actual filter code %
Y2=fft(detrend_rolling_vector);
Y2=filter_fft_vector(Y2); % frequency domain filtering in .oct file
y2=real(ifft(Y2));
% end of filter code %

smoothed_output(ii,1)=y2(L)+intercept+(L-1)*slope;

endfor

I have utilised compiled .oct files where possible to speed up the function.

[Jaroslav Hajek-2]

On Tue, Oct 20, 2009 at 2:58 PM, babelproofreader
<babelproofreader@...> wrote:

>
>>Have you tried implementing this in Octave code (.m)?
>
> Yes, I have; the code is below.
>
> function smoothed_output=fastrollingfftlowpass(data,fre,wid)
> %function y=fftlowpass(rolling_vector,fre,wid)
> if(exist('wid')~=1) wid=0.2; endif
> if(exist('fre')~=1) fre=0.05; endif
> n=rows(data);
> smoothed_output=data;
> L=16;
> L_2=L/2;
> rolling_vector=zeros(L,1);
> detrend_rolling_vector=zeros(L,1);
>
> for (ii=L:n)
>
> rolling_vector(1:L,1)=data(ii-(L-1):ii,1);
>
> % detrend the rolling_vector %
> intercept=rolling_vector(1,1);
> slope=(rolling_vector(L,1)-intercept)/(L-1);
> detrend_rolling_vector=detrend_vector(rolling_vector); % calls an .oct file
> "detrend_vector"
> % end of detrend code %
>
> % the actual filter code %
> Y2=fft(detrend_rolling_vector);
> Y2=filter_fft_vector(Y2); % frequency domain filtering in .oct file
> y2=real(ifft(Y2));
> % end of filter code %
>
> smoothed_output(ii,1)=y2(L)+intercept+(L-1)*slope;
>
> endfor
>
> I have utilised compiled .oct files where possible to speed up the function.
>
>
>
>

Instead of doing it in a for loop, I'd suggest you construct a matrix
with one rolling vector per column:

rmat = horzcat (cellslices (data, 1:n-L+1, L:n){:});

then, calculate your detrended vectors (this requires your function
can work column-wise):
dmat = detrend_matrix (rmat);
inter = rmat(1,:);
slope = (rmat(L,:) - inter) / L-1;

and do the filtering, utilizing the fast column-wise operation mode of fft:
y2 = fft (dmat);
y2=filter_fft_matrix (y2); # again you need to support column-wise operation
y2 = real (ifft (y2));

smoothed_output = (y2(L,:) + inter + (L-1)*slope)(:);

this will likely give you a much better performance.
If data is very long, you can avoid wasting memory by performing the
above operation on chunks of data,
chunk = 32000;
for i = L:chunk:N
  range = i-L+1:min(N,i+chunk-1);
  ## use the above operation with data(range) and smoothed_output(range)
  ## ...
endfor

hth

--
RNDr. Jaroslav Hajek
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[babelproofreader]

>Instead of doing it in a for loop, I'd suggest you construct a matrix
>with one rolling vector per column:

How would this compare, in terms of performance, to my function as it currently stands if I were able to successfully compile the whole function into a single .oct function using loops? I ask this question as I envisage that teaching myself how to/and rewriting/changing the .oct files that are currently called to operate over matrix columns will be a considerable effort.

[Jaroslav Hajek-2]

On Wed, Oct 21, 2009 at 1:25 PM, babelproofreader
<babelproofreader@...> wrote:

>
>>Instead of doing it in a for loop, I'd suggest you construct a matrix
>>with one rolling vector per column:
>
> How would this compare, in terms of performance, to my function as it
> currently stands if I were able to successfully compile the whole function
> into a single .oct function using loops? I ask this question as I envisage
> that teaching myself how to/and rewriting/changing the .oct files that are
> currently called to operate over matrix columns will be a considerable
> effort.
>
>

You won't really know until (if) you try.
I would expect my approach to perform better than the code you've
shown, but that's just my guess.
Also, if you're stuck with the old 3.0.1 version, things may be very
different from 3.3.50+ / 3.2.x.

hth

--
RNDr. Jaroslav Hajek
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[babelproofreader]

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.  

[David Grundberg-2]

babelproofreader skrev:

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

I had a quick look and I can't find anything obvious.  Also, it modifies
the value - the output is not the same as the input.  There's a lot of
intermediate results.  Any of these could be wrong.  Have you tried
printing out the result in each stage and comparing it to the correct
value from the m-code?

David
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[babelproofreader]

>...it modifies the value - the output is not the same as the input.

How so? When I plot the input against the output the two plots are identical.

>Have you tried printing out the result in each stage and comparing it to the correct
value from the m-code?

I compiled the code in stages and so I am sure that it is the filter code that is the problem. The code without the filtering, i.e.

#include <octave/oct.h>
#include <octave/dColVector.h>
#include <octave/parse.h>

DEFUN_DLD (pass3, 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 ();
  ifft_result = feval ("ifft", octave_value (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;
}

works OK.

[babelproofreader]

>I compiled the code in stages and so I am sure that it is the filter code that is the problem

The problem seems to be that the loop over the complex column vector output of fft fails. The simple code below

#include <octave/oct.h>
#include <octave/dColVector.h>
#include <octave/parse.h>
#include <math.h>

DEFUN_DLD (justchange, args, , "Help String")
{
  octave_value retval;
  ComplexColumnVector justchange_input = args(0).complex_column_vector_value ();

   for (octave_idx_type ii (0); ii < args(0).length (); ii++)
       {
       justchange_input(ii) = justchange_input(ii) * 1 ;
       }

  retval = justchange_input;

  return retval;
}

fails to compile, with error message

justchange.cc:13: error: no match for ‘operator=’ in ‘justchange_input.ComplexColumnVector::<anonymous>.MArray<std::complex<double> >::<anonymous>.Array<T>::operator() [with T = std::complex<double>](ii) = operator*(const octave_value&, const octave_value&)(((const octave_value&)(& octave_value(1))))’
/usr/include/c++/4.3/complex:1224: note: candidates are: std::complex<double>& std::complex<double>::operator=(double)
/usr/include/c++/4.3/complex:1159: note:                 std::complex<double>& std::complex<double>::operator=(const std::complex<double>&)

and this is essentially what I am trying to do.