[OctDev] pwelch.m bug report and patch

[<rol...@ie...>]

                                Hello,
I hope this is the right place to report bugs in octave-forge.

pwelch.m does not handle complex row vectors correctly.  This is
demonstrated in the script below.  White noise is filtered by a filter
with complex impulse response.  The filter frequency response peaks at
f/fs = 0.2 and so does the spectrum of the output signal.  However, if
the signal is a row vector, pwelch and tfe yield a peak at f/fs = 0.8.
The patch below (for pwelch.m as of octave-forge 2001.11.02, as found
in Debian octave-forge   2001.11.02-4) corrects this.  I believe that
the problem exists also in the code as of a few days ago, but I cannot
easily test it because I run rather old versions of Octave.
                                        hope this helps
                                               Roland

P.S. I am not (yet) subscribed to octave-dev, please CC me. Thanks.

----------- script to demonstrate the problem -----------
## pwelch_bug.m
## Created: 2004-05-11T10:48:25+0200  rol...@ie...
##   This code is put in the public domain.
## Demonstrate bug (incorrect handling of complex row vectors) in
## pwelch.m

Nresp = 60;                     # impulse response length
Wresp = 0.1;                    # filter bandwidth
h = exp(-Wresp * (0:(Nresp-1)));
Nfft = 256;
[H, w] = freqz(h, 1, Nfft);     # Frequency response

## Shift the filter response to F0.
F0 = 0.2;                       # centre frequency
hs = h .* exp(2i * pi * F0 * (0:(Nresp-1)));
Hs = freqz(hs, 1, Nfft);        # Frequency response

multiplot(2, 2);

## Generate input signal
Nsig = 10000;                   # signal length
x = 0.5 * randn(1, Nsig);
y = filter(hs, 1, x);

[Syy, Faxis]  = pwelch(y, Nfft, 1, 'whole'); # PSD of y
[Syyc, Faxis] = pwelch(y .', Nfft, 1, 'whole'); # PSD of y

Sxy  = tfe(x, y, Nfft, 1, 'whole');
Sxyc = tfe(x, y .', Nfft, 1, 'whole');

subplot(221);
clearplot;
title('Filter Magnitude');
xlabel('f/fs');
ylabel('20 log|H(f)| [dB]');
plot(w/pi, 20*log10(abs(Hs)), ';Filter response;');
hold on
plot(Faxis, 10*log10(abs(Syyc)), ';Syy col;');
plot(Faxis, 20*log10(abs(Sxyc)), '+;Sxy col;');

#subplot(222);
#clearplot;
#title('Filter Magnitude, row');
#xlabel('f/fs');
#ylabel('20 log|H(f)| [dB]');
#plot(w/pi, 20*log10(abs(Hs)), ';Filter response;');
plot(Faxis, 10*log10(abs(Syy)), ';Syy row;');
plot(Faxis, 20*log10(abs(Sxy)), '+;Sxy row;');

## Check filter phase
subplot(222);
clearplot;
title('Filter Phase, column');
xlabel('f/fs');
ylabel('arg(H(f)) [rad]');
plot(Faxis, unwrap(arg(Hs)), ';Filter phase;');

## Note that the cross-spectral density Sxy and the
## filter frequency response H are related as
## Sxy(f) = conj(H(f)):
## 
## Sxy(f) = Sxx(f) H^*(f)
## Syy(f) = Sxy(f) H(f) = Sxx(f) |H(f)|^2
## 
## Athanasios Papoulis
## Signal Analysis
## isbn 0-07-048460-0
## McGraw-Hill, 1984
## 9.3 Spectral Analysis, p. 312

plot(Faxis, -unwrap(arg(Sxyc)), '+;tfe phase;');
## plot(Faxis, -unwrap(arg(Sxy)), ';tfe phase row;');

subplot(223);
clearplot;
title('Estimate error, column');
xlabel('f/fs');
ylabel('20 log|H_est(f)/H(f)| [dB]');
plot(w/pi, 20*log10(abs(Sxyc .')) - 20*log10(abs(Hs)), ';Magnitude error;');

subplot(224);
clearplot;
title('Estimate error, column');
xlabel('f/fs');
ylabel('arg(H_est(f)/H(f)) [rad]');
plot(w/pi, -unwrap(arg(Sxyc .'))-unwrap(arg(Hs)), ';Phase error;');

## pwelch_bug.m ends here.


--------- patch against pwelch.m as of octave-forge 2001.11.02 ---------
*** pwelch.m	Wed Jan 23 03:05:25 2002
--- pwelch.m	Thu May 13 11:44:08 2004
***************
*** 99,104 ****
      error ([calledby, " data must be a vector"]);
    end
!   if columns(x) != 1, x = x'; end
!   if columns(y) != 1, y = y'; end
    if !isempty(y) && rows(x)!=rows(y)
      error ([calledby, " x and y vectors must be the same length"]);
--- 99,105 ----
      error ([calledby, " data must be a vector"]);
    end
!   ## Use simple transposition to avoid conjugation of complex x, y
!   if columns(x) != 1, x = x .'; end
!   if columns(y) != 1, y = y .'; end
    if !isempty(y) && rows(x)!=rows(y)
      error ([calledby, " x and y vectors must be the same length"]);

----------------- patch ends above this line ----------------------