pkg load signal clear all; N=20; fs=44100; f=1000; w=2*pi*f; t=0:1/fs:1 % 1 second % RAW square wave w_sqr = square(w*t); % RAW Square wave figure; plot(t,w_sqr); axis([0 1 -1.2 1.2]); grid on; hold on; audiowrite('test_square_wave_RAW.wav',w_sqr,fs); % approximation for t=0:1/fs:1 for n=1:N %Fourier-Fejér taper sum = sum + ((n-2*N-2)*(cos(pi*n)-1)*sin(n*w*t))/(2*pi*n*(N + 1)); end Ffej(i)=2*sum; i=i+1; sum=0; end; Ffej = Ffej'; plot(t,Ffej); axis([0 1 -1.2 1.2]); audiowrite('test_square_wave_Ffej.wav',Ffej,fs); figure; NFFT = 4096; L=length(w_saw); X = fft(w_saw,NFFT); X = X(1:NFFT/2+1); Px1 = X.*conj(X)/(NFFT*L); f = fs*(0:NFFT/2)/NFFT; plot(f,10*log10(Px1),'g', 'linewidth', 3); hold on; L=length(Ffej); X = fft(Ffej,NFFT); X = X(1:NFFT/2+1); Px2 = X.*conj(X)/(NFFT*L); f = fs*(0:NFFT/2)/NFFT; plot(f,10*log10(Px2),'k','linestyle','--', 'linewidth', 2); title('Single Sided Power Spectral Density'); xlabel('Frequency (Hz)') ylabel('Power Spectral Density- P_{xx} dB/Hz'); legend('w_sqr', 'Ffej'); ylim([-60 0]); grid on;to approximate square wave is really sloooow... .

I've tried to vectorize code but, (as being a old-school-programmer) can't figure out the vectorization logic ... just get errors. Is the sum calculation vectorizable and how to do it correctly?

N=20; fs=44100; f=1000; w=2*pi*f; t=0:1/fs:1 % approximation %Fourier-FejÃ©r taper sum_ = sum((((1:N)-2*N-2).*(cos(pi*(1:N))-1).*sin((1:N).*w.*(0:1/fs:1)'))./(2*pi*(1:N)*(N + 1)), 2); Ffej = 2*sum_; Ffej = Ffej';

Thank you very much!

My for-loop code took over one hour to get result (also seemed to use one core only) ... your code gave result immediately.

You could probably simplify/speedup the line by doing proper matrix algebra instead of just dot-products and sums, but that would require spending 5 minutes to actually understand what the code does. I'd rather not :-)

Just curious why Octave or Matlab tool makers can't automate that job and convert user slow loops to vectors. I have seen some field matlab designs full of slow loops that took hours to finish if it didn't break down.

@kaz I assume it's because Matlab is an interpreted language, so there isn't a compiler doing that kind of optimization.

Well I mean just a utility to convert loops to arrays. Surely they can do that just like they have done so many utilities with sub-licences and get more money as well...