% [freq CSM] = makeCSM(samples, fSample, nFFT)
%
% Bob Dougherty 4/8/2010
% computes the Cross Spectral Matrix. Similar to CSM_IQ_FFT in OptiNav beamforming software,
% but about 30 times slower.
%
% inputs:
% samples is a 2D array of sample values. The number rows is the number of times, and the
% number of columns is the number of channels (nMics)
% fSample is the sampling rate
% nFFT is the transform length. Perhaps 1024 or 2048. Ususally a power of 2.
%
% outputs:
% freq is an array of nFFT/2 - 1 analysis frequencies. freq(j) = j*(fSample/nFFT).
% CSM is the nMics x nMics x nFreq complex cross spectral matrix. It is normalized like a spectrum,
% not like a cross spectral density. If you want the cross spectral density normalization, divide by
% the analysis bandidth: fSample/nFFT.
%
function [freq CSM] = makeCSM(samples, fSample, nFFT)
nMics = size(samples,2);
nFreq = nFFT/2 - 1;
df = fSample/nFFT;
freq = zeros(nFreq);
for iFreq=1:nFreq
freq(iFreq) = iFreq*df;
end
CSM = zeros(nMics,nMics,nFreq);
for m=1:nMics
for n = m:nMics
Pxy = cpsd(samples(:,n),samples(:,m),hanning(nFFT),nFFT/2,nFFT,fSample);
for iFreq=1:nFreq
CSM(n,m,iFreq) = Pxy(iFreq+1)*df;
if (m~=n)
CSM(m,n,iFreq) = CSM(n,m,iFreq)';
end
end
end
end