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...
