FFT and changing frequency and vectorizing for loop

classic Classic list List threaded Threaded
5 messages Options
RT
Reply | Threaded
Open this post in threaded view
|

FFT and changing frequency and vectorizing for loop

RT
FFT and changing frequency and vectorizing for loop

Greetings All

I can increase and decrease the frequency of the
signal using the combination of fft and a series expansion for loop in the code below
but if the signal/array is to large it becomes extremely
slow (an array that's 1x44100 takes about 2 mins to complete) I'm sure it has to do with the for loop but
I'm not exactly sure how to vectorize them to improve performance

Any recommendations

tia sal22 

%create signal
clear all, clc,clf,tic
x= linspace(0,2*pi,44100)'; 

%Used in exporting to ycalc audio file make sure in sync with above
freq_orig=1;
freq_new=4
vertoff=0;
vertoffConj=0;
vertoffInv=0;
vertoffInvConj=0;
phaseshift=(0)*pi/180 ; %can use mod to limit to 180 degrees

y=sin(freq_orig*(x)); 
[size_r,size_c]=size(y);

N=size_r; %to test make 50
T=2*pi;
dt=T/N;
t=linspace(0,T-dt,N)';
phase = 0;
f0 = 1/T; % Exactly, one period 

y=(y/max(abs(y))*.8)/2; %make the max amplitude here
C = fft(y)/N; % No semicolon to display output
 

A = real(C); 
B = imag(C)*-1; %I needed to multiply by -1 to get the correct sign

% Single-Sided (f >= 0) 
An = [A(1); 2*A(2:round(N/2)); A(round(N/2)+1)]; %needed to put the ' to get vaules in rows 
Bn = [B(1); 2*B(2:round(N/2)); B(round(N/2)+1)]; %needed to put the ' to get vaules in ro

pmax=N/2;
ycalc=zeros(N,1); %preallocating space for ycalc
w=0;
for p=2:pmax
                %
        %%1 step) re-create signal using equation 
         ycalc=ycalc+An(p)*cos(freq_new*(p-1).*t-phaseshift)+Bn(p)*sin(freq_new*(p-1).*t-phaseshift)+(vertoff/pmax);
         w=w+(360/(pmax-1)); %used to create phaseshift
         phaseshift=w;
        
end;
fprintf('\n- Completed in %4.4fsec or %4.4fmins\n',toc,toc/60);

subplot(2,1,1), plot(y),title('Orginal Signal'); 
subplot(2,1,2),plot(ycalc),title('FFT new signal'); 


--
-

_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: FFT and changing frequency and vectorizing for loop

Sergei Steshenko

--- On Sat, 4/2/11, Rick T <[hidden email]> wrote:

From: Rick T <[hidden email]>
Subject: FFT and changing frequency and vectorizing for loop
To: [hidden email]
Date: Saturday, April 2, 2011, 3:02 PM

FFT and changing frequency and vectorizing for loop
Greetings All
I can increase and decrease the frequency of thesignal using the combination of fft and a series expansion for loop in the code below

but if the signal/array is to large it becomes extremelyslow (an array that's 1x44100 takes about 2 mins to complete) I'm sure it has to do with the for loop butI'm not exactly sure how to vectorize them to improve performance


Any recommendations
tia sal22 
%create signalclear all, clc,clf,ticx= linspace(0,2*pi,44100)'; 
%Used in exporting to ycalc audio file make sure in sync with above

freq_orig=1;freq_new=4vertoff=0;vertoffConj=0;vertoffInv=0;vertoffInvConj=0;phaseshift=(0)*pi/180 ; %can use mod to limit to 180 degrees


y=sin(freq_orig*(x)); [size_r,size_c]=size(y);
N=size_r; %to test make 50T=2*pi;dt=T/N;t=linspace(0,T-dt,N)';phase = 0;

f0 = 1/T; % Exactly, one period 
y=(y/max(abs(y))*.8)/2; %make the max amplitude hereC = fft(y)/N; % No semicolon to display output 
A = real(C); 

B = imag(C)*-1; %I needed to multiply by -1 to get the correct sign
% Single-Sided (f >= 0) An = [A(1); 2*A(2:round(N/2)); A(round(N/2)+1)]; %needed to put the ' to get vaules in rows 

Bn = [B(1); 2*B(2:round(N/2)); B(round(N/2)+1)]; %needed to put the ' to get vaules in ro
pmax=N/2;ycalc=zeros(N,1); %preallocating space for ycalcw=0;for p=2:pmax

                %        %%1 step) re-create signal using equation          ycalc=ycalc+An(p)*cos(freq_new*(p-1).*t-phaseshift)+Bn(p)*sin(freq_new*(p-1).*t-phaseshift)+(vertoff/pmax);

         w=w+(360/(pmax-1)); %used to create phaseshift         phaseshift=w;        end;fprintf('\n- Completed in %4.4fsec or %4.4fmins\n',toc,toc/60);


subplot(2,1,1), plot(y),title('Orginal Signal'); subplot(2,1,2),plot(ycalc),title('FFT new signal'); 

--
-


-----Inline Attachment Follows-----

_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave

I don't think you are doing a right thing - not WRT loop/vectorization,
but WRT final result.

If I understand correctly your code, you multiply by 4 each individual
frequency. If so, your signal will appear in your output 4 times in a row.

Try a simple example with, say, just 16 samples.

There is a package for 'octave' doing polyphase resampling:

http://octave.sourceforge.net/signal/function/resample.html
.

Regards,
  Sergei.

_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
RT
Reply | Threaded
Open this post in threaded view
|

Re: FFT and changing frequency and vectorizing for loop

RT
It works it's just slow due to the for loop here's an example of the output which is correct


On Sat, Apr 2, 2011 at 2:41 PM, Sergei Steshenko <[hidden email]> wrote:

--- On Sat, 4/2/11, Rick T <[hidden email]> wrote:

From: Rick T <[hidden email]>
Subject: FFT and changing frequency and vectorizing for loop
To: [hidden email]
Date: Saturday, April 2, 2011, 3:02 PM

FFT and changing frequency and vectorizing for loop
Greetings All
I can increase and decrease the frequency of thesignal using the combination of fft and a series expansion for loop in the code below

but if the signal/array is to large it becomes extremelyslow (an array that's 1x44100 takes about 2 mins to complete) I'm sure it has to do with the for loop butI'm not exactly sure how to vectorize them to improve performance


Any recommendations
tia sal22 
%create signalclear all, clc,clf,ticx= linspace(0,2*pi,44100)'; 
%Used in exporting to ycalc audio file make sure in sync with above

freq_orig=1;freq_new=4vertoff=0;vertoffConj=0;vertoffInv=0;vertoffInvConj=0;phaseshift=(0)*pi/180 ; %can use mod to limit to 180 degrees


y=sin(freq_orig*(x)); [size_r,size_c]=size(y);
N=size_r; %to test make 50T=2*pi;dt=T/N;t=linspace(0,T-dt,N)';phase = 0;

f0 = 1/T; % Exactly, one period 
y=(y/max(abs(y))*.8)/2; %make the max amplitude hereC = fft(y)/N; % No semicolon to display output 
A = real(C); 

B = imag(C)*-1; %I needed to multiply by -1 to get the correct sign
% Single-Sided (f >= 0) An = [A(1); 2*A(2:round(N/2)); A(round(N/2)+1)]; %needed to put the ' to get vaules in rows 

Bn = [B(1); 2*B(2:round(N/2)); B(round(N/2)+1)]; %needed to put the ' to get vaules in ro
pmax=N/2;ycalc=zeros(N,1); %preallocating space for ycalcw=0;for p=2:pmax

                %        %%1 step) re-create signal using equation          ycalc=ycalc+An(p)*cos(freq_new*(p-1).*t-phaseshift)+Bn(p)*sin(freq_new*(p-1).*t-phaseshift)+(vertoff/pmax);

         w=w+(360/(pmax-1)); %used to create phaseshift         phaseshift=w;        end;fprintf('\n- Completed in %4.4fsec or %4.4fmins\n',toc,toc/60);


subplot(2,1,1), plot(y),title('Orginal Signal'); subplot(2,1,2),plot(ycalc),title('FFT new signal'); 

--
-


-----Inline Attachment Follows-----

_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave

I don't think you are doing a right thing - not WRT loop/vectorization,
but WRT final result.

If I understand correctly your code, you multiply by 4 each individual
frequency. If so, your signal will appear in your output 4 times in a row.

Try a simple example with, say, just 16 samples.

There is a package for 'octave' doing polyphase resampling:

http://octave.sourceforge.net/signal/function/resample.html
.

Regards,
 Sergei.




--


_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: FFT and changing frequency and vectorizing for loop

Sergei Steshenko
In reply to this post by RT

--- On Sat, 4/2/11, Rick T <[hidden email]> wrote:

From: Rick T <[hidden email]>
Subject: Re: FFT and changing frequency and vectorizing for loop
To: "Sergei Steshenko" <[hidden email]>
Cc: [hidden email]
Date: Saturday, April 2, 2011, 7:00 PM

It works it's just slow due to the for loop here's an example of the output which is correcthttp://dl.dropbox.com/u/6576402/questions/q1.png




On Sat, Apr 2, 2011 at 2:41 PM, Sergei Steshenko <[hidden email]> wrote:





--- On Sat, 4/2/11, Rick T <[hidden email]> wrote:



From: Rick T <[hidden email]>

Subject: FFT and changing frequency and vectorizing for loop

To: [hidden email]

Date: Saturday, April 2, 2011, 3:02 PM



FFT and changing frequency and vectorizing for loop

Greetings All

I can increase and decrease the frequency of thesignal using the combination of fft and a series expansion for loop in the code below



but if the signal/array is to large it becomes extremelyslow (an array that's 1x44100 takes about 2 mins to complete) I'm sure it has to do with the for loop butI'm not exactly sure how to vectorize them to improve performance








Any recommendations

tia sal22 

%create signalclear all, clc,clf,ticx= linspace(0,2*pi,44100)'; 

%Used in exporting to ycalc audio file make sure in sync with above



freq_orig=1;freq_new=4vertoff=0;vertoffConj=0;vertoffInv=0;vertoffInvConj=0;phaseshift=(0)*pi/180 ; %can use mod to limit to 180 degrees





y=sin(freq_orig*(x)); [size_r,size_c]=size(y);

N=size_r; %to test make 50T=2*pi;dt=T/N;t=linspace(0,T-dt,N)';phase = 0;



f0 = 1/T; % Exactly, one period 

y=(y/max(abs(y))*.8)/2; %make the max amplitude hereC = fft(y)/N; % No semicolon to display output 

A = real(C); 



B = imag(C)*-1; %I needed to multiply by -1 to get the correct sign

% Single-Sided (f >= 0) An = [A(1); 2*A(2:round(N/2)); A(round(N/2)+1)]; %needed to put the ' to get vaules in rows 



Bn = [B(1); 2*B(2:round(N/2)); B(round(N/2)+1)]; %needed to put the ' to get vaules in ro

pmax=N/2;ycalc=zeros(N,1); %preallocating space for ycalcw=0;for p=2:pmax



                %        %%1 step) re-create signal using equation          ycalc=ycalc+An(p)*cos(freq_new*(p-1).*t-phaseshift)+Bn(p)*sin(freq_new*(p-1).*t-phaseshift)+(vertoff/pmax);



         w=w+(360/(pmax-1)); %used to create phaseshift         phaseshift=w;        end;fprintf('\n- Completed in %4.4fsec or %4.4fmins\n',toc,toc/60);





subplot(2,1,1), plot(y),title('Orginal Signal'); subplot(2,1,2),plot(ycalc),title('FFT new signal'); 



--

-





-----Inline Attachment Follows-----



_______________________________________________

Help-octave mailing list

[hidden email]

https://mailman.cae.wisc.edu/listinfo/help-octave



I don't think you are doing a right thing - not WRT loop/vectorization,

but WRT final result.



If I understand correctly your code, you multiply by 4 each individual

frequency. If so, your signal will appear in your output 4 times in a row.



Try a simple example with, say, just 16 samples.



There is a package for 'octave' doing polyphase resampling:



http://octave.sourceforge.net/signal/function/resample.html

.



Regards,

  Sergei.






--

No, your example doesn't prove it works. I.e. the image shows one original
'sin' period, and four 'sin' periods in the "upsampled" signal.

Take as input signal, for example, a triangle pulse - you'll see
four triangle pulsed after your "upsampling". So, if you want to upsample
a song this way, you'll hear it four times.

What you did is simple replication in time domain - can be done quickly
using just several lines of code without any FFT in the first place.

Regards,
  Sergei.

P.S. Please do not use HTML;
P.P.S. Please do no top post - this list prefers bottom posting.
_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: FFT and changing frequency and vectorizing for loop

Pertti Paakkonen
In reply to this post by RT
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