DFT takes too long to run
1 view (last 30 days)
Show older comments
Anonymous
on 31 Oct 2021
Answered: Sulaymon Eshkabilov
on 31 Oct 2021
I am creating a PSD plot from a DFT after analyzing an audio file. It is taking so long to run. Is there a way to shorten the run time?
% Plotting Audio to Find Sample Rate
[y,Fs] = audioread('guitar.wav');
fprintf('Sampling Rate = %f samples/s', Fs) %sampling rate
samples=[1,length(y)-(3.5*Fs)]; %cut of the dissipated section of audio
[y1,Fs] = audioread('guitar.wav',samples);
format shortg;
dt = 1/Fs;
t = 0:dt:(length(y1)*dt)-dt; % time values
t2 = transpose(t);
figure(1)
plot(t,y1);
xlabel('Seconds');
ylabel('Amplitude');
T = max(t2); %sampling period, period (length of sampling in DFT) [s]
w = 2*pi/T; % angular frequency (based on *entire* sample time for DFT)
L = length(t2); % number of samples
% Discrete Fourier Transform (DFT)
a0 = 2/L*sum(y1); % constant term
N = floor(L/2); % maximum of floor(L/2): interesting/necessary part
a = zeros(N,1); b = a; % pre-allocate
Y = a0/2*ones(1,L); %allocate space for Y
for n = 1:N % find the first N coefficeients
a(n) = 2/L*sum( y1.*cos(w*n*t2) ); % sum over *samples*
b(n) = 2/L*sum( y1.*sin(w*n*t2) ); % sum over *samples*
Y = Y + a(n)*cos(n*w*t2) + b(n)*sin(n*w*t2); % estimated function
end
% Plot PSD
W = w*(1:N)'; % frequency range
P = sqrt(a.^2 + b.^2); % calculate power
figure(2); % Plot PSD in rad/s
plot(W,20*log10(P),'ko-'); % plot power [dB] against frequency [rad/s]
grid on; axis tight; set(gca,'fontsize', 14);
xlabel('\omega [rad/s]','fontsize', 14); ylabel('Power [dB]','fontsize', 14);
% Plot PSD in Hz
figure(3);
plot(w/(2*pi)*(1:N),20*log10(P),'ko-'); % plot power [dB] against frequency [Hz]
grid on; axis tight; set(gca,'fontsize', 14);
xlabel('f [Hz]','fontsize', 14); ylabel('Power [dB]','fontsize', 14);
0 Comments
Accepted Answer
Sulaymon Eshkabilov
on 31 Oct 2021
First of all it looks like that the heavy calc burden is occuring within [for .. end] loop. Therefore, a(n) and b(n) should be computed without a loop just using a vectorization. E.g.:
n = 1:N;
[NN, T2]=meshgrid(n,t2);
a = 2/L*sum(y1.*cos(w*NN.*T2));
b = 2/L*sum( y1.*sin(w*NN.*T2));
Y = sum(a.*cos(NN*w.*T2) + b.*sin(NN*w.*T2));
0 Comments
More Answers (0)
See Also
Categories
Find more on Parametric Spectral Estimation in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!