function [pkso] = Get_Audio_PSD(fname, intwid)
A = audioread(fname);
Ainfo = audioinfo(fname);
Fs = Ainfo.SampleRate;
disp(['Processing audio at ', num2str(Fs./1000), 'kHz'])
try
B = A(1E5:2E5,1);
C = A(1E5:2E5,2);
catch
B = A(:,1);
C = A(:,2);   
end
[ab,b] = xcorr(B);
[ac,b] = xcorr(C);

x = ab;
N = length(x);
xdft = fft(x);
xdft = xdft(1:N/2+1);
psdx = (1/(Fs*N)) * abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);
freq = 0:Fs/length(x):Fs/2;

% plot(freq,10*log10(psdx))

psdxo(:,1) = psdx;
xdo(:,1) = freq;
% hold on

x = ac;
N = length(x);
xdft = fft(x);
xdft = xdft(1:N/2+1);
psdx = (1/(Fs*N)) * abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);
freq = 0:Fs/length(x):Fs/2;

%plot(freq,10*log10(psdx), 'r')
% grid on
% title(fname)
% xlabel('Frequency (Hz)')
% ylabel('Power/Frequency (dB/Hz)')
% legend('Channel 1', 'Channel 2')

psdxo(:,2) = psdx;
xdo(:,2) = freq;

% From here we only deal with psdxo and xdo...

[~, xl] = max(psdxo);
disp(xl)
funda = (xdo(xl(1)));
disp(['Found fundamental at ', num2str(funda)])

%Integrate the peaks over the harmonics based on the funamental position...
size(psdxo)
% for f = 1:floor(length(xdo)./xl)
for f = 1:(floor(length(psdxo)./xl)-1)
    ceil(f*xl+intwid/2)
   pkso(f, 1) = 10*log10(max(psdxo(floor(f*xl-intwid/2):ceil(f*xl+intwid/2), 1)));
   pkso(f, 2) = 10*log10(max(psdxo(floor(f*xl-intwid/2):ceil(f*xl+intwid/2), 2)));
end

