Search code examples
matlaboctavesox

converting a vector loop to a for loop because of memory limitations


I have a piece of Octave / Matlab code that I received a lot of help from Andy and the group with. The issue i’m running into now is not enough memory to make signals that are longer in duration.

My plan to work around this is:

1) convert the vector loop into a for loop. (having a problem here)
2) have the for loop export each segment of the loop as a wav file instead of doing what the vector code does which is append it. (having a problem here)
3) join each wave file segment using sox.

Most of the examples online go from for loops to vectorized loops not the other way around, any ideas? I’m also open to other recommendations to fix my memory problem. Note: I’m working with 1 gig of ram on a rasberry pi 2 and it works and it’s pretty fast, I’m just trying to get signals with a longer duration, and exporting each segment should allow for this.

I’m using Octave which is compatible with Matlab.

See the working vectorized code below: Its based off of Paul Nasca stretch algorithm found here http://www.paulnasca.com/algorithms-created-by-me#TOC-PaulStretch-extreme-sound-stretching-algorithm

urlwrite('http://www.onewithall.net/rttmpfiles/3sec8000.wav','3sec8000.wav'); 
inputfn='3sec8000.wav' %change this to test another file
[d, fs, bps]=wavread(inputfn);
inputlen=rows  (d)/fs;

printf ("Original duration of file in seconds = %.2f s\n", rows (d)/fs);

dur=60; %duration / length you want the file to be in seconds
stretch = dur/rows(d)*fs; %how much I need to stretch the file to get it to be the duration I want
windowsize = round (0.25 * fs);

step = round ((windowsize/2)/stretch);

## original window
fwin = @(x) (1-x.^2).^1.25;
win = fwin (linspace (-1, 1, windowsize));

#win = hanning (windowsize)';

## build index
ind = (bsxfun (@plus, 1:windowsize, (0:step:(rows(d)-windowsize))'))';
cols_ind = columns(ind);

## Only use left channel
left_seg = d(:,1)(ind);
clear d ind;

## Apply window
left_seg = bsxfun (@times, left_seg, win');

## FFT
fft_left_seg = fft (left_seg);
clear left_seg

#keyboard

## overwrite phases with random phases
fft_rand_phase_left = fft_left_seg.*exp(i*2*pi*rand(size(fft_left_seg)));
clear fft_left_seg;

ifft_left = ifft (fft_rand_phase_left);
clear fft_rand_phase_left;

## window again
ifft_left = bsxfun (@times, real(ifft_left), win');

## restore the windowed segments with half windowsize shift
restore_step = floor(windowsize/2);
ind2 = (bsxfun (@plus, 1:windowsize, (0:restore_step:(restore_step*(cols_ind-1)))'))';
left_stretched = sparse (ind2(:), repmat(1:columns (ind2), rows(ind2), 1)(:), real(ifft_left(:)), ind2(end, end), cols_ind);
clear ind2 ifft_left win;

left_stretched = full (sum (left_stretched, 2));

## normalize
left_stretched = 0.8 * left_stretched./max(left_stretched);
printf ("converted %s =%.2f(s) file to stretched.wav = %.2f(s)\n", inputfn, inputlen, rows (left_stretched)/fs);
wavwrite (left_stretched, fs, bps, "streched.wav");

I tried tracking the problem down by putting in display('line') at key points. And it looks like the line

left_stretched = sparse (ind2(:), repmat(1:columns (ind2), rows(ind2), 1)(:), real(ifft_left(:)), ind2(end, end), cols_ind);

The above line seems to have a problem only when I run out of ram. It says error subscript indices must be either positive integers or logicals. please note this only happens when I run out of memory when trying to use long durations by setting dur=60*1800. If I set dur=60*10 everything works.


Solution

  • Do you remember me? I'm the author of the initial code which you've posted. Below the code as for loop. I've tested this with an output len of 800s.

    ## based on http://hypermammut.sourceforge.net/paulstretch/
    ## https://github.com/paulnasca/paulstretch_python/blob/master/paulstretch_steps.png
    
    more off
    inputfn = "original.wav"
    [d, fs, bps] = wavread (inputfn);
    inputlen=rows  (d)/fs;
    
    printf ("Original  duration of file in seconds = %.2f s\n", rows (d)/fs);
    
    target_duration = 60; # in seconds
    stretch = target_duration/inputlen;
    # 1/4 s window len
    windowsize = round (0.25 * fs);
    
    # stepwidth between windows
    step = round ((windowsize/2)/stretch);
    numsteps = floor((rows(d)-windowsize)/step);
    
    ## restore the windowed segments with half windowsize shift
    restore_step = floor (windowsize / 2);
    
    ## stetched duration
    stretched_len = (numsteps*restore_step+windowsize)/fs;
    printf ("Stretched duration of file in seconds = %.2f s\n", stretched_len);
    stretched = zeros (numsteps*restore_step+windowsize, 1);
    if (!exist ("out", "dir"))
      mkdir ("out");
    endif
    
    ## Matrix which holds the freq of the maximum amplitude and the max. amplitude
    chunks_stats = zeros (numsteps, 2);
    
    ## original window
    fwin = @(x) (1-x.^2).^1.25;
    win = fwin (linspace (-1, 1, windowsize));
    
    ## loop over all windows
    for k = 1:numsteps
      if (! mod(k, 50))
        printf ("Calculating chunk %i of %i...\n", k, numsteps);
        fflush (stdout);
      endif
    
      ## Only use left channel
      s_ind = (k - 1) * step + 1;
      e_ind = s_ind + windowsize - 1;
      tmp = win' .* d(s_ind:e_ind, 1);
    
      ## FFT, overwrite phases with random phases and IFFT
      tmp = fft(tmp);
      [m, ind] = max (abs(tmp(1:numel(tmp)/2)));
      # Frequency in Hz
      chunks_stats(k, 1) = (ind-1)/windowsize*fs;
      # max Amplitude
      chunks_stats(k, 2) = m;
      printf ("Freq =  %.2f Hz, max = %.2f\n", chunks_stats(k, :));
      tmp = ifft (tmp .* exp(i*2*pi*rand(size(tmp))));
    
      ## window again
      tmp = win' .* real (tmp);
      fn = sprintf ("out/out_%04i.wav", k);
      wavwrite (tmp, fs, bps, fn);
      s_ind = (k - 1) * restore_step + 1;
      e_ind = s_ind + windowsize - 1;
      stretched (s_ind:e_ind) += tmp;
    endfor
    
    ## normalize
    stretched = 0.8 * stretched./max(stretched);
    wavwrite (stretched, fs, bps, "stretched.wav");
    

    If you want to write multiple wavs to concatenate them later it's a little bit more difficult because the overlapping windows. But I think this code will run fine on a BeagleBoneBlack.

    EDIT: Added saving chunks to separate files and collect maximum amplitude and frequency of this signal per chunk to chunk_stats.