DFT takes too long to run

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

Answer 1

[Sulaymon Eshkabilov]

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