Possible conflict with the Signal Processing Toolbox in MATLAB.

[vikaasa]

I recently observed a mismatch in output, where two machines that ran the same code utilizing this public domain package, had a significant difference in output (but no warning or error message).

The reason was one of the machines had the Signal Processing Toolbox installed, which was listed as a dependency on running the matlab.codetools.requiredfilesandproducts function.

We had observed that the machine with the Signal Processing Toolbox installed and listed as a dependency, gave an incorrect output. On the other hand, the machine with the correct output had only base MATLAB as a dependency.

[nicolas-cellier-aka-nice]

Hi, I think that this is more or less intentional: if you take a look at https://github.com/metocean/diwasp/blob/ae31ec24b36f77ad435fbac5da42f44604434a84/private/diwasp_csd.m on line 10 you will see a test for existence of cpsd function. This is a function of the signal toolbox which computes cross-spectrum of two signals. So if you have the signal toolbox, you take this if branch, else you take the else one.

The thing is that the two branches are not equivalent... In the cpsd branch, there is no windowing function (we pass nfft as the window parameter which just means that we want a rectangular Hamming window of length nfft). In the else branch, there is a Hann window See https://en.wikipedia.org/wiki/Hann_function.

Let's choose some meaningless x,y signals:

fs=40; t=(1:512)'/fs;
x=cos(11*t+0.07*t.^2).*sin(2.5*t)+0.2*randn(size(t));
y=cos(11*t+1+0.05*t.^2).*cos(2.4*t-1)+0.2*randn(size(t));
nfft=64;

then compute cross spectrum with the DIWASP way (else branch):

my_hann=0.5*(1-cos(2*pi*(1:nfft/2)/(nfft+1))); %use my_ann so as to not overload hann function
win = [my_hann my_hann(end:-1:1)];
nw = length(win);
nseg=fix(length(x)/nw);
S = zeros(nfft,1);
for iseg=0:nseg-1
    ind=nw*iseg+[1:nw];
    xw = win'.*x(ind);
    yw = win'.*y(ind);
    Px = fft(xw,nfft);
    Py = fft(yw,nfft);
    Pxy = Py.*conj(Px);
    S = S + Pxy;
end
nfac=(fs*nseg*norm(win)^2);
S=[S(1); 2*S(2:nfft/2); S(nfft/2+1)]/nfac;
f=(fs/nfft)*[0:nfft/2]';

versus the cpsd way:

[S1, f1] = cpsd(x,y,nfft,0,nfft,fs);

or modify invocation of cpsd to use the same window as DIWASP:

[S2, f2] = cpsd(x,y,win,0,nfft,fs);

then plot the magnitude of the 3 cross spectrum above:

figure; plot(f,abs(S),'r.',f1,abs(S1),'b*',f2,abs(S2),'ro');
legend('DIWASP','cpsd original','cpsd modified');
saveas(gcf,'diwasp_csd','png');

we can see that the original cpsd (blue stars) differs on the peak, due to the window, but the modified one matches DIWASP (red circles and dots)...

So one thing that could be changed is to compute the Hann window outside the if, then use it in both branches.

Personnally, I prefer to use an overlap like is the default of cpsd:

noverlap=nfft/2;
[S1, f1] = cpsd(x,y,nfft,noverlap,nfft,fs);

[nicolas-cellier-aka-nice]

Wait, there is more:

Not only the magnitude is different in the two branches, but also the angle!

figure; plot(f,angle(S),'r.',f1,angle(S1),'b*',f2,angle(S2),'ro');
legend('DIWASP','cpsd original','cpsd modified');
saveas(gcf,'diwasp_csd_angle','png');

You can see in obtained figure that angles of cpsd branch (circles and stars) are exact opposites of the DIWASP branch (dots)...

Why? because the convention in cpsd(x,y) is usingfft(X) .* conj(fft(y)) while the DIWASP branch uses the conjugate: Py*conj(Px)...

This clearly qualifies as a bug IMO.