bugGNU Octave - Bugs: bug #50820, [octave forge] (signal) pwelch...

 
 

bug #50820: [octave forge] (signal) pwelch incorrectly discards DC component

Submitter:  - <_102371>
Submitted:  Mon 17 Apr 2017 09:39:47 PM UTC
   
 
Category:  Octave Package Severity:  1 - Wish
Priority:  1 - Later Item Group:  Matlab Compatibility
Status:  Patch Submitted Assigned to:  None
Originator Name:  Open/Closed:  * Open
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Mon 04 Mar 2019 08:19:22 PM UTC, comment #26: 

I'm attaching another update to pwelch, this time so that if a second output is requested but fs is not supplied, the second output will contain a normalized frequency vector.  That seems to be what is documented for Matlab and the previous version for Octave, but for me, the Octave version failed in this case.



(file #46418)

John W. Eaton <jwe>
Group administrator
Wed 09 Jan 2019 10:30:35 PM UTC, comment #25: 

OK.  Is there anything I can do to help that happen?

John W. Eaton <jwe>
Group administrator
Wed 09 Jan 2019 09:52:13 PM UTC, comment #24: 

I want to at least get a signal package 1.4.1 out with some smaller bug fixes that have accumulated since 1.4.0, and save this rewrite for 2.0.0.

Mike Miller <mtmiller>
Group Member
Wed 09 Jan 2019 09:46:06 PM UTC, comment #23: 

I updated the docstring.  My latest version is attached.

Mike: would it be OK to include this version in the signal package now?  I think this version is far more compatible with pwelch in the current Matlab signal processing toolbox.

(file #45919)

John W. Eaton <jwe>
Group administrator
Thu 11 Oct 2018 06:09:58 PM UTC, comment #22: 

I think I've integrated your latest changes in my version.  Please check the attached file.

If you have more changes, please make them relative to this version.  Otherwise, I think it will be unnecessarily difficult to merge your changes given all the other modifications I've made.

Yes, full name and email.

(file #45183)

John W. Eaton <jwe>
Group administrator
Wed 10 Oct 2018 12:55:48 AM UTC, comment #21: 

Sorry for the quick posts, forgot to mention that MATLAB's default for "short" vector is strictly greater than 8. Any shorter and it throws an error. I will leave it to you to decide if Octave should enforce the same as well. I am personally indifferent as I do not see many cases where such a short vector will be used.

- <_102371>
Tue 09 Oct 2018 11:52:30 PM UTC, comment #20: 

P.P.S enbw() is a function in the signal processing toolbox in MATLAB. We can add that single line of code to the Octave signals toolbox as well. Here is an example function.

(file #45174)

- <_102371>
Tue 09 Oct 2018 11:26:39 PM UTC, comment #19: 

You are most certainly welcome :). What is the usual credit format? Name? Email? Both? Other?

P.S. If there are other things in the signal processing side of things that need help, do let me know. I could certainly lend my signal processing expertise and contribute where I am able to.

- <_102371>
Tue 09 Oct 2018 11:23:14 PM UTC, comment #18: 

John: Further testing revealed a couple of bugs that I fixed. Here is my file. Please merge your changes to this file, if that is not too much work for you. If not, let me know and I will do that when I find time from my research.

Fixed:

  • Trace arguments interaction with confidence estimation. Apparently regardless of the trace arguments, confidence values remain the same. So the average periodogram needs to be saved regardless of the trace type, which I have done so.
  • Confidence interval is returned for n_ffts = 1, which is also fixed.
  • "mean" is a valid trace argument. Fixed that in parse_args()



(file #45173)

- <_102371>
Tue 09 Oct 2018 07:14:52 PM UTC, comment #17: 

Thanks for the reply.  I made some style changes and switched to using the second arg to fft as well as indexing + reshape to handle N-dimensional arrays.  My current version is attached.

Please let me know how you'd like to be credited in the commit message.

Thanks again for your help.


(file #45171)

John W. Eaton <jwe>
Group administrator
Tue 09 Oct 2018 03:36:39 PM UTC, comment #16: 

John: I noticed that as well, but for the time being I did not focus on optimising existing code; just get all the features to work. If you are happy with my changes to the code, pwelch now is completely compatible with MATLAB, barring few tests.

1) Yes, fft() pads zero as required and we really do not require that before hand.
2) I was going to go for reshape as well initially, it was really a flip of a coin for this one. I just needed something that worked.
3) I cannot think of any case where it would make sense to subtract 3. Speaking of which, we should include a sanity check to make sure the specified size or window vector does not exceed the signal size.
4) I am currently not sure how MATLAB enforces short vector size, I have to test it out next time I am in the university. But on top of my head, I think that there should at least be 3 elements for shifting/onesided PSD to work properly.

- <_102371>
Tue 09 Oct 2018 03:10:54 PM UTC, comment #15: 

Thanks for the update.

Doesn't fft(X,N) pad X with zeros as needed?  It seems easier to do that than to insert X into a zero-padded XX array.  That's why I changed the call to use fft(X,N).

Instead of generating a cell array of colons, you can use indexing like this:


x = reshape (1:60, [10,2,3])  ## 3d array for illustration.
xx = x(start:end,:)  ## all columns in all pages as a 2d array.
reshape (xx, [end-start+1,2,3])  ## IF you need 3d array again.


I was also confused by the subtraction of 3.  It was in the original pwelch.m from Octave Forge.  Are there any cases where it is needed to avoid generating a window size that is too large?

What should happen for short vectors?  Is there a minimum size limit that should be enforced?

John W. Eaton <jwe>
Group administrator
Mon 08 Oct 2018 08:43:10 PM UTC, comment #14: 

John: Here is what I have so far. I have NOT used any code from MATLAB. The only MATLAB resource I have used is their documentation, to see what arguments the functions use and their definitions.
I have commented the changes with my initials and made relevant references where necessary. They can be removed if you want to.

The only issue so far is that: confidence interval computation requires functions from statistics package. I have currently made it so that it errors out if the package is not installed/not loaded. I do not know what should be done here, the function is not simple enough to reimplement without too much headaches.

- <_102371>
Mon 08 Oct 2018 05:29:47 PM UTC, comment #13: 

STPR: could you post the changes you have made?  Also, it's important that any changes you make are your work and not derived in any way from Matlab code that you may have access to.


John W. Eaton <jwe>
Group administrator
Sat 06 Oct 2018 03:14:45 AM UTC, comment #12: 

Here is an update on this issue. I tested the function with some EMG data I had used extensively in my past research and found some big and small issues (matrices not supported, wrong indices shifted, some arguments {psd_type, trace arguments, etc.}).

I decided to fix the issues to get things going; and managed to fix several. I am currently finalising some changes with the last argument (confidence args). Once I change that and test it to make sure it gives the same outputs as MATLAB does, I will post the code here.

- <_102371>
Thu 04 Oct 2018 01:28:24 PM UTC, comment #11: 

I can certainly test it, I am part of the signal processing group at my university. I have a bunch of signals I can test out to see if I observe something odd.
I also see a fixme section. Do you need descriptions for that? I can certainly provide them if needed.

- <_102371>
Wed 03 Oct 2018 06:21:32 PM UTC, comment #10: 

I'm attaching an updated version of pwelch.  I wouldn't be surprised to find that I've missed something, but as far as I can tell at the moment, it is compatible with the current Matlab pwelch behavior.

There were MANY features in the previous version that no longer work.  I actually think that's a good thing because this function is complicated enough as it is.

Some testing is definitely needed.  I could use help with that from someone who actually knows something about signal processing and uses pwelch.

(file #45133)

John W. Eaton <jwe>
Group administrator
Mon 01 Oct 2018 03:25:58 PM UTC, comment #9: 

It looks like spectrogram also has a similar set of arguments.

John W. Eaton <jwe>
Group administrator
Fri 28 Sep 2018 11:27:57 PM UTC, comment #8: 

Thanks for starting on this. Yeah the argument handling looks like more than half the work of this function. Looks good to me so far.

I just noticed that the Matlab function 'periodogram' has significant overlap in its command line arguments with 'pwelch'. Maybe once we have something working there could be some reuse.

Mike Miller <mtmiller>
Group Member
Fri 28 Sep 2018 12:20:58 PM UTC, comment #7: 

OK, I started working on this.  Just doing the option parsing and setting the default option and parameter values, and validating the options and parameters for the current pwelch behavior is quite complicated bu itself without including any extra features.  I think I have most of it now.  My current version is attached.

Note that making it strictly compatible with the current pwelch interface in Matlab will drop quite a few features from the current version in Octave's signal package.  I'm personally OK with that because it will reduce the maintenance burden.  But users may not like that their scripts suddenly fail with an upgrade to a new signal package version.

Comments and suggestions welcome.


(file #45112)

John W. Eaton <jwe>
Group administrator
Tue 25 Sep 2018 01:51:25 AM UTC, comment #6: 

I would vote for #1. I really think the aim should be to be compatible with the latest behaviour on MATLAB, which I think should be the goal of Octave. Most people I know who use MATLAB almost all use close to the latest version (most universities have licenses to use the new ones anyway). And as Mike pointed out, people who want the old interface can download the older version.

- <_102371>
Fri 21 Sep 2018 06:24:46 PM UTC, comment #5: 

With comment #3, I was in favor of your option 1, throw out all the complications and aim for compatibility with the latest Matlab behavior.

I don't think there's a way to keep the old behavior without breaking the interface in some way anyway. If someone really wants the old pwelch, I would rather point them to an old release or the hg history.

Mike Miller <mtmiller>
Group Member
Fri 21 Sep 2018 04:40:24 PM UTC, comment #4: 

I would like to work on updating pwelch to match the current interface and calculations done by the same function in Matlab's signal processing toolbox.

What would be your preference for backward compatibility?  I see a couple of possibilities:

1. Throw out all the complications of handling old Matlab behavior and target the latest interface.

2. Same as 1, but also keep the existing function as old_pwelch (or something).  If someone really wants the old behavior, they can copy this function file to a directory in their personal path.

3. Keep the old function and automatically forward to it if someone has set the internal "compatib" variable.  Otherwise, just try to match current behavior.  This approach seems friendlier at minimal cost.  I would further mark the old code as unmaintained and make no promise for attempting to update it.  It seems to me that the new code will be easier to maintain if we don't attempt to mix in the old in the same function file.

Comments and suggestions?  I can begin working on the compatible interface and argument processing code and identify any computations that need to be added or updated to match current behavior.

John W. Eaton <jwe>
Group administrator
Wed 21 Mar 2018 01:57:14 AM UTC, comment #3: 

Octave's pwelch already supports updated Matlab defaults, but via a strange "R11" or "R12" argument passed in.

I haven't really modified pwelch myself, but it might make sense to try to simplify all the code that tries to be compatible with different versions of Matlab and just go with the latest. It's hard enough keeping up with the changes that Matlab makes, but trying to keep code optionally backwards compatible with older versions doesn't seem worth it.

So someone might want to look at the "compatib" internal variable in pwelch.m and the various interpretations of positional arguments that is already done, and choose whichever is closest to the current Matlab way.

Mike Miller <mtmiller>
Group Member
Mon 17 Apr 2017 10:10:05 PM UTC, comment #2: 

Thanks for reply.
I suspected as much, it looked as if the DC was being subtracted. Adding "no-strip" argument results in the expected figure.

Yes, MATLAB does indeed not remove DC from the input. As a suggestion, would it be possible to remove the DC if a "strip" option is specified rather than keep if "no-strip" is specified? It will make things more uniform across Octave and MATLAB.

- <_102371>
Mon 17 Apr 2017 09:59:41 PM UTC, comment #1: 

I think the issue here is that in octave's pwelch the default
is to subtract mean from the signal, and in matlab's apparently not.

Try to add "no-strip" option:

[X,f]=pwelch(x,fs,0,fs,fs, "no-strip");

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Mon 17 Apr 2017 09:39:47 PM UTC, original submission:  

Current setup:
Windows 10, GNU Octave 4.2.1, Signal Package 1.3.2.

Procedure to replicate issue:

Generate any signal (or use existing one if desired) with a non-zero DC. Use the pwelch() method using one second window, no overlap to get the power spectrum. Observe the DC coefficient on the result.

What to expect?:

The DC term to be the actual average of the signal.

What happens?:

The DC term is zero.

Additional info:

I used the attached code to test for the issue. I generated a sinusoid with some DC offset. As a comparison, I generated the same graph on MATLAB. MATLAB has the correct DC term while Octave does not. See the attached graphs to observe the results. The graph name *_ML.png was generated in MATLAB. This shows the proper DC term. Octave does not.

- <_102371>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #46418:  pwelch.m added by jwe (20KiB - text/x-objcsrc)
file #45919:  pwelch.m-2019-01-09 added by jwe (20KiB - application/octet-stream)
file #45183:  pwelch.m-2018-10-11 added by jwe (14KiB - application/octet-stream)
file #45174:  enbw.m added by _102371 (1019B - application/octet-stream)
file #45173:  pwelch.m added by _102371 (17KiB - application/octet-stream)
file #45171:  pwelch.m-2018-10-09 added by jwe (13KiB - application/octet-stream)
file #45165:  pwelch.m added by _102371 (15KiB - text/x-objcsrc)
file #45133:  pwelch.m-2018-10-03 added by jwe (12KiB - application/octet-stream)
file #45112:  pwelch.m added by jwe (7KiB - text/x-objcsrc)
file #40424:  pwtest.m added by _102371 (259B - application/octet-stream)
file #40425:  PW_Test_ML.png added by _102371 (23KiB - image/png)
file #40426:  PW_Test_Oct.png added by _102371 (14KiB - image/png)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by jwe (Updated the item)
  • -email is unavailable- added by dasergatskov (Posted a comment)
  • -email is unavailable- added by _102371 (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 22 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2019-03-04 jwe Attached File- Added pwelch.m, #46418
    2019-02-26 mtmiller Severity2 - Minor 1 - Wish
        StatusPostponed Patch Submitted
        Releaseother dev
    2019-01-09 jwe Attached File- Added pwelch.m-2019-01-09, #45919
    2018-10-11 jwe Attached File- Added pwelch.m-2018-10-11, #45183
    2018-10-09 _102371 Attached File- Added enbw.m, #45174
    2018-10-09 _102371 Attached File- Added pwelch.m, #45173
    2018-10-09 jwe Attached File- Added pwelch.m-2018-10-09, #45171
    2018-10-08 _102371 Attached File- Added pwelch.m, #45165
    2018-10-03 jwe Attached File- Added pwelch.m-2018-10-03, #45133
    2018-09-28 jwe Attached File- Added pwelch.m, #45112
    2018-09-26 mtmiller Carbon-CopyRemoved 80942 -
    2018-03-21 mtmiller Severity3 - Normal 2 - Minor
        Priority5 - Normal 1 - Later
        Item GroupIncorrect Result Matlab Compatibility
        StatusNone Postponed
        Release4.2.1 other
    2017-08-12 jwe Summarypwelch incorrectly discards DC component [octave forge] (signal) pwelch incorrectly discards DC component
    2017-04-17 _102371 Attached File- Added pwtest.m, #40424
        Attached File- Added PW_Test_ML.png, #40425
        Attached File- Added PW_Test_Oct.png, #40426

    Back to the top

    Powered by Savane 3.13-4448.
    Corresponding source code