Frequency response of FM modulation/demodulation chain with phase derivative demodulation

Question

Frequency modulating a carrier by white noise and then demodulating the complex signal using discrete derivative of phase it appears that the discriminator is acting as a low-pass filter.

How do I calculate the rolloff vs modulating frequency for a given phase-based discriminator's impulse response?

I know about Octave's freqz function, but applying it to the discriminator's impulse response gives the frequency response of the discriminator, which in my understanding would translate to amplitude distortion of the demodulated signal, so this is not directly the frequency response of the modulate->demodulate chain which I'm looking for.

pkg load signal;
orig_mono_240Khz = randn(1,150000);
orig_mono_240Khz = orig_mono_240Khz - mean(orig_mono_240Khz);
max_val = max([abs(max(orig_mono_240Khz)), abs(min(orig_mono_240Khz))])*1.2;

subplot(2,1,1);
plot(abs(fftshift(fft(orig_mono_240Khz))));
axis("tight");

phase_changes = j*75000*2*pi.*orig_mono_240Khz/max_val;
signal = resample(exp(cumsum(phase_changes/240000)), 10, 1);
signal = signal - min(signal);
signal = round(signal/max([max(imag(signal)),max(real(signal))])*255);


signal = signal - mean(signal);
rcv_240 = decimate(signal, 10);
unwrp_phase = unwrap(angle(rcv_240));
unwrp_phase = unwrp_phase - mean(unwrp_phase);
disc_resp = [1, 0, -1];
phase_drv = conv(unwrp_phase, disc_resp)(100:end-100);
phase_drv = phase_drv-mean(phase_drv);

subplot(2,1,2);
plot(abs(fftshift(fft(phase_drv))));
axis("tight");

Answer 1

You are essentially seeing a Cascade-Integrator-Comb (CIC) response which is identical to a moving average filter (Aliased Sinc function magnitude response) as seen with CIC filter structures.

Consider what is happening in units of phase:

You start with a white noise signal which is translated from magnitude directly to units of frequency in the FM modulation process.

To demodulate you convert to phase (which is an integration of the instantaneous frequency of your waveform, since phase is the integral of frequency) so if you evaluated the frequency response of the resulting phase versus time signal you would expect to see the low-pass frequency response of an integrator, rolling off at 1/f. (You aren't actually doing any integration, and extracting the phase itself from your waveform is not filtering, the phase is already filtered as if you had integrated your original modulation waveform: For example, the phase for white-FM rolls off at 1/f in frequency). Also be aware that for constant envelope waveforms (no AM) the power spectral density is directly proportional to the phase fluctuations-- for small angles, the power of the sidebands is given directly by the phase versus time.

You then compute the differences in phase between samples over N samples to recover the modulating signal-- this is a digital approximation of differentiation and is the classic "Comb Filter", as with N large the magnitude response represents the teeth of a comb.

$$1-z^{-N}$$

So essentially the process matches the cascade of an accumulator (integrator) with a comb, which is identical to a moving average filter over N samples.

To see how these are equivalent:

Moving average as a geometric series: $$ \begin{align} & 1+z^{-1}+z^{-2}+z^{-3} + ... + z^{-N-1} \\ & = \sum_{n=0}^{N-1}z^{-n} \\ & = \frac{1-z^{-N}}{1-z^{-1}} \end{align}$$

accumulator: $\frac{1}{1-z^{-1}}$

comb: $1-z^{-N}$

Also interesting is the poles and zeros of each which also explain how the moving average is identical to a CIC: The comb filter has N total zeros starting at z = 1 and equally spaced on the unit circle (the roots of $1-z^{-N}$). The integrator has a pole at z=1. A moving average filter has N-1 zeros, all at the same locations as the comb filter except for z=1.

From your spectrum it looks like you are using N=1. As you increase N you would expect the envelope of the result to approach a Sinc function consistent with a moving average filter.