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)

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.

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)

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)

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)

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.

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.

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.

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)

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.

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.

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.

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.

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.