Trying to do DFT and getting 'Unable to perform assignment because the left and right sides have a different number of elements' error.

I'm trying create a PSD plot of a recorded guitar string pluck (wav file), manually, by using a DFT. However, I keep getting the error 'Unable to perform assignment because the left and right sides have a different number of elements' when the code reaches 'a(n)'. I'm not sure how to fix this error. Any help is super appreciated.

[y,Fs] = audioread('guitar.wav');
fprintf('Sampling Rate = %f samples/s', Fs) %sampling rate

format shortg;
y(1:2.5*Fs,:) = []; % this deletes the matrix elements corresponding to the index range 1:2.5*fs
                    % acts as a way to 'zoom in' to pure tone
y(end-1*Fs:end,:) = [];
dt = 1/Fs;
t = 0:dt:(length(y)*dt)-dt; % time values

figure(1)
plot(t,y);
xlabel('Seconds');
ylabel('Amplitude');

T = max(t);   %sampling period, period (length of sampling in DFT) [s]
w = 2*pi/T;     % angular frequency (based on *entire* sample time for DFT)
L = length(t);          % number of samples

% Discrete Fourier Transform (DFT)
a0 = 2/L*sum(y);        % 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( y.*cos(w*n*t) );            % sum over *samples*
    b(n) = 2/L*sum( y.*sin(w*n*t) );            % sum over *samples*
    Y = Y + a(n)*cos(n*w*t) + b(n)*sin(n*w*t);  % 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

[Srijith Kasaragod]

The expression "2/L*sum( y.*cos(w*n*t) )" returns an array of size 1xM. Equating it to "a(n)" is bound to produce an error as "a" is an Nx1 vector. So the LHS of the equation expects a single element but is being passed a 1xM array. You can resolve the issue by modifying 'zeros' function in assigning "a" as:

a=zeros(N,M); 

Regards.