function [xdo, psdxo] = Get_Audio_PSD(fname, forcefreq)
Ainfo = audioinfo(fname);
Fs = Ainfo.SampleRate;

if nargin==2
    Fs = forcefreq;
end
A = audioread(fname);

disp(Fs)
B = A(:,1);
C = A(:,2);  


[ab,b] = xcorr(B);
[ac,b] = xcorr(C);

x = ab;
x = x((length(x)./3):(length(x)./4)*2);
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);
disp(Fs)
disp(N)
freq = 0:(Fs/length(x)):(Fs/2);



plot(freq,10*log10(psdx))
disp(['freq', num2str(size(freq))])
disp(['psdx', num2str(size(psdx))])


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

x = ac;
x = x((length(x)./3):(length(x)./4)*2);

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;


