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

[Jaroslav Hajek]

On Tue, Oct 20, 2009 at 2:58 PM, babelproofreader
<address@hidden> 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