FFT y cambio de frecuencia y vectorización del bucle FOR

Question

Puedo aumentar y disminuir la frecuencia de una señal utilizando la combinación de fft y un bucle FOR de expansión de la serie de Fourier en el código de abajo, pero si la señal/array es muy grande se vuelve extremadamente lento (un array que es 1x44100 tarda unos 2 minutos en completarse) Estoy seguro de que tiene que ver con el bucle for pero no estoy exactamente seguro de cómo vectorizarlo para mejorar el rendimiento. Hay que tener en cuenta que esto se utilizará con señales de audio de 3 a 6 mins. El array 1x44100 es sólo un segundo y tarda unos 2 mins en completarse

Alguna recomendación

%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)]; 
Bn = [B(1); 2*B(2:round(N/2)); B(round(N/2)+1)];

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');

Aquí hay una foto de la trama si alguien quiere ver la salida, que es correcta el bucle FOR es sólo muy, muy lento

Answer 1

Creo que puedes sustituir la parte principal del bucle for por la siguiente línea

ycalc=An.*cos(freq_new*(1:(p-1))*t-phaseshift) +Bn.*sin(freq_new*(1:(p-1))*t-phaseshift)+(vertoff/pmax);