TransferFunction.FrequencyResponse return half sided spectrum + 1 (Nyquist + 1)

[ToGoOrNotToGo]

While review and learning the code from the TransferFunction class, I wonder why FrequencyResponse returns a half-sided spectrum + 1 (Nyquist + 1), but, for example, GroupDelay not. Is there a reason or is it a bug?

FrequencyResponse:
real.Take(length / 2 + 1)

GroupDelay:
var gd = new double[fftSize / 2];

[ToGoOrNotToGo]

Oh and the for index in GroupDelay should start at 0, I think...
for (var i = 1; i <= gd.Length; i++)
->
for (var i = 0; i < gd.Length; i++)

otherwise you will get a peak in the last index.

[ar1st0crat]

No, it's not a bug; it's just how I've implemented GD. The reason for fftSize / 2 + 1, I suppose, is quite obvious: due to DFT properties, this number of first coefficients carry all information - the rest are symmetric. For example, if the spectrum is 0 1 2 3 4 3 2 1, then 3 2 1 can be discarded in order to save memory, and we keep only the first 8/2 + 1 = 5 samples. That's why PowerSpectrum / MagnitudeSpectrum / FrequencyResponse return exactly fftSize/2 + 1 samples.

Group delay, by definition, is the derivative (see Wikipedia, for example). It can be naively implemented like here:

group_delay = -diff(unwrap(angle(h))) / diff(w)

diff function in Python returns derivative of size n - 1 where n is the size of array:

a = [5, 6, 7, 8, 9]
diff(a)

# returns [1, 1, 1, 1]

Essentially, we're just excluding zero frequency element from computations, that's why the length is smaller by 1, i.e.fftSize / 2. Like in the example above, the first resulting value is not 5 (a[0] - 0), but 1 (a[1] - a[0]).

PS. In fact, the implementation is slightly more difficult - see here, but the 'derivative thing' is the same.

[ToGoOrNotToGo]

Yes. The half-sided spectrum is clear. But why + 1? Look/debug into the result for e.g. Bessel:
Visual Studio QuickWatch:

	Name	Value	Type
	[252]	-0.79615712452371667	double
	[253]	-0.7959939873750681	double
	[254]	-0.79589693674270023	double
	[255]	4	double

For me the last number (4) looks strange...

[ar1st0crat]

I've already explained why +1 (see my previous answer). Seems like you're asking about the group delay that has size fftSize / 2 in current implementation (without +1). Yeah, I see the problem with the last sample, but it has nothing to do with indexing. Indexing is OK and the algorithm is implemented properly. I've checked other filters with different parameters and sometimes these edge effects are gone, sometimes they happen at the beginning. Right now I don't know why it's happening. Hope I'll figure it out...

[ar1st0crat]

Well, according to Rick Lyons:

Recently I was using this group delay estimation method on a 6th-order IIR bandpass filter and my results appeared to be incorrect. As digital filter guru Sanjeev Sarpal pointed out to me, the above group delay estimation method works correctly for FIR filters but not always correctly for IIR filters. In fact, MATLAB's 'grpdelay()' command uses the above scheme for FIR filters, but uses a scheme called "the Shpak algorithm" for IIR filters. So, to be safe, don't use the above scheme for IIR filters.

So for better results for IIR filters I'll check out the Shpak algorithm...

[ToGoOrNotToGo]

OK thank you. I understand now fftSize / 2 + 1 in FrequencyResponse.
I was confused because the different handling in group delay. And because the LinePlot didnt display the last (strange ;-) ) peak from group delay, so maybe there is an issue in LinePlot too!?

[ar1st0crat]

Quite possible )) LinePlot is not a part of NWaves, so I don't really care about it ). I wrote it very quickly just to have general visual idea about what's goin' on

[ar1st0crat]

I've written a complete analog of sciPy group_delay function (and maybe the MATLAB grpdelay as well - I didn't test it) with default values for w and whole parameters. It still has problems (the algorithm is essentially the same, i.e. it's not the Shpak algorithm), but in general works better and at least it's unified with existing implementations in standard libraries. Closin' this, since I think for now it's more than enough, and I'll be working on other features.