Re: FFT and changing frequency and vectorizing for loop

[Pertti Paakkonen]

Rick,

Probably vectorization does not help here. This is how your for p=2:pmax
loop is vectorized (with comments on some oddities):

p = (2:pmax)';
%Fix An and Bn to correct size. The first and last elements in the original
code were never used.
An = 2*A(p); 
Bn = 2*B(p); 
phaseshift = 360/(pmax-1)*(p-2); %Why not in radians?
M  = (An*ones(1,N)).*cos(freq_new*(p-1)*t'-phaseshift*ones(1,N));
M += (Bn*ones(1,N)).*sin(freq_new*(p-1)*t'-phaseshift*ones(1,N)); %Matlab
incompatible notation
ycalc2 = sum(M)+(pmax-1)*vertoff/pmax; %Take vertoff/pmax out from
summation.

Result is the same with maximum error of 5e-15.

Vectorization involves expanding your frequency components and time vector
to at least five 22049 times 44100 matrices with 7.2 GB each. Unless you
have a huge physical memory this will fail. Even with smaller time vectors
the looping is more efficient than vectorization. This is probably a memory
handling issue because number of sin() and cos() executions remain the same
and with very small time vectors (say, less than 1000) vectorized code runs
faster.

In my test Rick's original looped code took 5.72 seconds to run while code
above took 8.16 seconds when N was "only" 8192. Case N=16328 ran out of
memory. System was octave 2.1.57-7 (i686-pc-linux-gnu) on a 12 GB quad-core.

If one finds better way to vectorze this please let us know :-)

Pertti


--
View this message in context: 
http://octave.1599824.n4.nabble.com/FFT-and-changing-frequency-and-vectorizing-for-loop-tp3422784p3423091.html
Sent from the Octave - General mailing list archive at Nabble.com.