Search code examples
matlabsignal-processingsimilarity

Comparison of Higher Order Ambisonics Signals

I'm working on a project with matlab where I have to compare two (Higher Order) Ambisonics impulse responses. I've calculated the FFT and plotted the smoothed energy of both signals over frequency:

Frequency plot

They seem to be correlated but I want to find a value to estimate the similarity. I can't use the function xcorr because my signals are Nx25 instead of Nx1.

Assume that a and b are my signals, both are 8000x25 matrices:

fs = 48e3;
fft1 = fft(a);
fft2 = fft(b);
l1 = length(a);
l2 = length(b);

% calculating frequencies 
f = fs*(0:(l1/2))/l1;

% calculating single sided spectrum for both signals

P1 = 10*log10(abs(fft1).^2);
P1 = P1(1:l1/2+1);
P1(2:end-1) = 2*P1(2:end-1);

P2 = 10*log10(abs(fft2).^2);
P2 = P2(1:l2/2+1);
P2(2:end-1) = 2*P2(2:end-1);

% plotting both spectra
figure; plot(f,smooth(P1,150)); hold on; plot(f,smooth(P2,150)); grid on

So I have 25 values per time (let's say each one of them represent a channel of my signal). So xcorr will allow me to compare one channel with another one. What I need to do is to compare the signal as a whole entity because every channel has spatial information that is important to the signal. So comparing one channel with the other won't lead me to an exact statement about the similarity.

I would appreciate your help and advice on how to compare these signals.

Solution

  • There are many ways to compare signals. Let's distinguish two main options first: (1) the signals should be compared at each time point, without shifting, or (2) the signals should be compared in a time-invariant manner, allowing a shift of the one with respect to the other.

    Case (2) typically comes down to a cross-correlation function, such as computed by xcorr. The value of the maximum of this function indicates the correlation under a time shift that maximizes the correlation. That is, we shift one signal in time to find the position where it maximally correlates with the other, then quantify that correlation. It can be computed as follows:

    sim = max(ifft(fft(a).*conj(fft(b))));

    The above computes the correlation for each channel of a and b independently, producing a value for each channel. Since you want to compare all channels of a to all channels of b, you need to set up a matrix. A simple way to do this is using the new implicit singleton expansion in MATLAB:

    b = reshape(b,[],1,25);
sim = max(ifft(fft(a).*conj(fft(b))));
sim = squeeze(sim);

    The above produces a 25x25 matrix of correlations. If you get an error computing it, you probably have an older version of MATLAB. This more ugly bit of code does the same thing:

    b = reshape(b,[],1,25);
sim = max(ifft(bsxfun(@times,fft(a),conj(fft(b)))));
sim = squeeze(sim);

    Case (1) is more interesting because there are so many more options besides correlation: (root) mean square error, mean absolute error, rank correlation, mutual information, total correlation, and a long list of etceteras. Some of these are trivial to compute, some are more complicated. But even the most complicated one you can compute for all combinations of channels using two for loops:

    sim = zeros(25,25);
for ii = 1:25
  for jj = 1:25
    sim(ii,jj) = similarity_function(a(:,ii),b(:,jj)); % whatever you want to fill in here!
  end
end

    It depends very much on the characteristic of the signals, which measure is going to be most useful. If you're interested in something simple like the root mean square error, it can be computed for each pair again using implicit singleton expansion:

    b = reshape(b,[],1,25);
sim = sqrt(mean((a-b).^2));
sim = squeeze(sim);

    Finally, we need to reduce the 25x25 matrix of pairwise channel comparisons to a single measure. Again, there are different options. One possibility is to simply look at the two most similar channels, or two most dissimilar channels:

    max(sim(:))
min(sim(:))

    or we can take some statistic related to all combinations:

    mean(sim(:))
median(sim(:))
sum(sim(:))
sqrt(mean(sim(:).^2));

    The other option is to find the channel in b that best matches each channel in a:

    sim = min(sim,[],2);

    This returns a value for each of the channels of a. Again, you can take the maximum or minimum value, or compute any other statistic as discussed earlier.