WFM Stereo Broadcasting Standards - OIRT vs CCIR

[radiolistener]

I’m currently trying to decode stereo from the test sample shared above. It seems I overlooked something in my earlier analysis. In CCIR, the DSB signal is transmitted without a carrier, so once you recover the carrier (which appears not so easy with proper phase), its makes L-R recovery straightforward. However, in OIRT the signal includes a carrier. When the carrier is locked with a PLL and down-converted, a noticeable DC offset appears. I’d prefer not to rely on a HPF filter. What would be the option to remove this offset/carrier?

A second complication is that this makes DSB amplitude recovery problematic, since it requires knowing the exact carrier level.

[Njk]

The subcarrier level is about 100 mVpp. It can be measured when silence is broadcasting

[radiolistener]

I was not referring to the absolute subcarrier amplitude. For proper stereo decoding it’s essential that the L+R and L–R paths maintain equal gain. The L+R component is straightforward, but the L–R signal first needs to be separated from the carrier and demodulated from DSB. This can be done with a synchronous detector using the PLL-locked carrier. However, because the carrier is present within DSB signal bandwidth, the result includes a DC offset and require HPF to remove it.

I’d prefer not to use a high-pass filter just to remove that DC. My idea was to try subtracting the carrier directly, by taking the PLL output and attenuating it by 14 dB, then substract it from composite. The issue is that the actual level of the carrier may vary between transmitters, so a fixed 14 dB attenuation would likely leave a residual offset due to this mismatch.

Any idea on how to extract L-R from DSB which include carrier but with no use of HPF?

[radiolistener]

I came across an interesting book by С.В.Мелихов, "Аналоговое и Цифровое Радиовещание" (S.V. Melikhov, "Analog and Digital Broadcasting") – see p.98 – which describes details of the MPX stereo composite for both OIRT polar modulation with a 31.25 kHz subcarrier and the CCIR system with a 19 kHz pilot tone. The book includes block diagrams for different demodulator types along with brief explanations.

It turns out that the OIRT composite is essentially equivalent to the CCIR composite, with one key difference: instead of a pilot tone, the carrier is present in the DSB signal, attenuated by 14 dB. In this form, stereo can indeed be decoded with just two diodes.

Here’s how a test stereo signal looks in Audacity (test composite is attached in flac file)

As you can see, the upper half-wave contains the left channel and the lower half-wave the right channel. I didn’t do anything special here – this effect occurs naturally from combining the L+R component with the DSB-C modulated by L–R component.

Below is an example block diagram of a simple polar stereo detector:

The description notes that the main advantage of the polar detector is its simplicity and ease of adjustment. The drawbacks are limited channel separation and increased distortion during detection, especially for higher audio frequencies. This limited separation occurs because the phase angle corresponding to zero crossings of the subcarrier (2θ′) is not constant at 180° for each half-period, but depends on the relative amplitudes of signals A and B.

The instability of the cutoff angles of diodes VD1 and VD2 leads to nonlinear distortion of signals A and B during detection. The THD of a polar detector in the low and mid-frequency range is about 2–3%. When detecting higher audio frequencies, the distortion increases further, reaching 6–8%. This is due to the fact that the subcarrier frequency of the polar-modulated signal is only 31.25 kHz / 15 kHz = 2.08 times higher than the upper audio frequency, meaning that within one period of the highest audio envelope there are only slightly more than two subcarrier cycles.

Because of these drawbacks, polar detectors are not used in high-quality receivers.

The book also presents a more refined sum-difference stereo decoder:

According to the description, thanks to full-wave detection of the supersonic portion of the polar-modulated signal – where more than four subcarrier periods fit into one period of the highest audio frequency envelope – nonlinear distortion does not exceed 1–2%. Channel separation depends on circuit balancing (adjusted by potentiometers R7 and R8) and can reach about 32–36 dB.

Another interesting design is a keyed stereo decoder that uses time-multiplexing:

It is noted that the shorter the keying pulses, the greater the channel separation and the more accurately the A and B signals are reproduced. To achieve >30 dB separation and <1% distortion over 0.3–5 kHz, the control pulse width must not exceed 5 µs.

Update: test-oirt-composite-50us.flac attachmend was updated due to resample aliasing issue in previous version

[radiolistener]

According to the book, the OIRT stereo composite consists of 0.8*(L+R) + AM-modulated 0.8*(L–R), with the AM carrier suppressed by 14 dB using a partial carrier suppression circuit (Q = 100 ± 5):


This can be expressed in Octave as:

Code: [Select]

composite = 0.8*add_LR + (0.8*sub_LR .* osc3125) + 0.2*osc3125;

I generated a test signal with these parameters (attached). Note that in my version I simply added the DSB-SC with a -14 dB (0.2) residual carrier at 31.25 kHz, whereas the book specifies that the transmitter uses an AM modulator and then suppresses the carrier by 14 dB. In theory, this results in some attenuation near DC, meaning low frequencies may be slightly reduced. In my test signal, that attenuation is absent, so bass components may sound louder. I think this is the meaning of K(f) function mentioned in the standard:

2.1.2.2 A signal S is produced equal to one half of the difference between signals A and B mentioned above. This signal, S, is pre-emphasized in the same way as signal M. The pre-emphasized signal, S, is used for the amplitude modulation of a sub-carrier at 31.25 kHz; the spectrum of the amplitude-modulated sub-carrier is formed so that the sub-carrier amplitude is reduced by 14 dB and the spectral components of the given modulating signal appear to be transformed as follows:

\$\overline{K}(f) = \frac{1 + j\,6.4 f}{5 + j\,6.4 f}\$

where f is equal to each frequency component (kHz).

The test file with an IQ recording containing an FM-modulated test composite (50 kHz deviation, 50 µs pre-emphasis) is attached. If anyone has the chance to test this signal on a receiver with an OIRT stereo decoder, I’d be grateful for feedback – in particular: are the left and right channels reproduced correctly, and is the loudness and frequency response as expected?

The file is provided in FLAC format to save space. FLAC is a lossless format, so no information is lost during conversion. Since most receivers only accept WAV, you can easily convert it using the sox utility by specifying the FLAC file as input and WAV as output:

Code: [Select]

sox test-oirt-wfm-50us.flac test-oirt-wfm-50us.wav

PS: The test file is generated with the following pipeline:
1) Input wav file has 48 kHz sample rate
2) Apply FIR LPF with cutoff=15 kHz and 140 dB stop-band attenuation (296 taps)
3) Upsample from 48 kHz to 192 kHz
4) Apply pre-emphasis approximated IIR: \$H(z) = 1 + k - k z^{-1};\$ where \$k = F_s \cdot \tau\$
5) composite = 0.8*add_LR + ((0.8*sub_LR) .* osc3125) + 0.2*osc3125;
6) Upsample from 192 kHz to 960 kHz
7) Do FM modulation with trapezoid integration:

Code: [Select]

    wfm_df = 50000;         % FM deviation Δf = ±50 kHz
    phi = 2*pi*wfm_df/Fs * cumtrapz(composite);     % Trapezoid integration,   H(s)=1/s -> H(z)=T/2 * (1+z^-1)/(1-z^-1)
    iq = exp( 1j*phi );
8) Apply FIR LPF with cutoff=115.5 kHz (50+31.25+15 kHz bandwidth + 19.25 kHz safeguard) and 140 dB stop-band attenuation (384 taps)
9) Downsample 960 kHz to 384 kHz
10) Save to flac with 8-bit resolution

Update: test-oirt-wfm-50us.flac was updated due to resample images issue in previous version

[radiolistener]

The book states that a keyed stereo decoder with time-multiplexing of the stereo channels can achieve better performance than a sum-difference stereo decoder. As I understand it, this implies that in DSP one could simply apply upsampling for more accurate phasing, then sample at the peaks of the oscillator locked by the PLL, and finally apply a low-pass filter to the positive and negative half-waves, resulting in a high-quality signal. Am I correct in this interpretation? How much should the sampling rate be increased in practice to achieve good channel separation and low THD+N?

I’m currently working on decoding using a PLL to recover the carrier, then mixing the signal down to the PLL-locked frequency, applying an LPF, integrating the DC component introduced by the subcarrier with 1-st order IIR LPF with 5 Hz cut-off, and subtracting it from the demodulated L–R. The L+R component is obtained simply by applying the same LPF to composite. The decoder works, but channel separation is somewhat not excelent, and a parasitic buzzing tone appears in the background around 2.5 kHz. A decoded file is attached.

I haven’t been able to determine the source of this tone. I noticed it seems related to the input sampling rate: at 192 kHz, the buzz is present, but when using Fs = 31.25 × 6 = 187.5 kHz, it disappears. It seems likely that PLL noise due to imperfect phase alignment with the samples is causing it. I’m using a PLL with an arctangent phase detector. What could be the cause? Interestingly, the same PLL works without noise when decoding a CCIR signal at Fs=192 kHz.

I experimented with adding a low-pass filter at the PLL phase detector input. This does shift the “humming” component to a lower frequency—roughly corresponding to the loop filter cutoff—but the noise still persists. Moreover, lowering the cutoff below 100 Hz leads to self-oscillation. What methods can be used to further reduce this background noise?

Update: A minor issue with the resampling process was discovered, which was introducing image components into the L+R band. I’ve fixed it and re-uploaded the composite, WFM, and decoded files. However, this did not affect the background noise during decoding — it is still present.

[Njk]

The test file with an IQ recording containing an FM-modulated test composite (50 kHz deviation, 50 µs pre-emphasis) is attached. If anyone has the chance to test this signal on a receiver with an OIRT stereo decoder, I’d be grateful for feedback – in particular: are the left and right channels reproduced correctly, and is the loudness and frequency response as expected?

Could you provide an example for such device?

[Njk]

BTW, the presence of subcarrier in the OIRT baseband signal provides one more advantage. In case of multi-path reception, everything becomes distorted, including the subcarrier sine wave. That facilitates easy indication of the multi-path condition. Use a tank tuned to the second harmonic of the subcarrier (62.5 kHz) to measure the harmonic's amplitude and show the number to the user. This is implemented in the L. RP-015 decoder (see the schematic diagram above).

[gf]

I’d prefer not to use a high-pass filter just to remove that DC. My idea was to try subtracting the carrier directly, by taking the PLL output and attenuating it by 14 dB, then substract it from composite. The issue is that the actual level of the carrier may vary between transmitters, so a fixed 14 dB attenuation would likely leave a residual offset due to this mismatch.

Any idea on how to extract L-R from DSB which include carrier but with no use of HPF?

Ideally, if the composite signal is exactly k1*(L+R)/2 + k2*subcarrier*(L-R)/2 + 0.19953 * subcarrier, where subcarrier has a peak amplitude of 1, then the DC offset (after down-conversion by multiplication with exp(-2i*pi*f_subcarrier*t)) should be a fixed constant value of 0.19953/2 = 0.099763, which you could subtract. However, in the case of imperfections, a high pass filter is likely the only option for removing an unknown DC component. You could use a very low cutoff, say 1 Hz. Note, by the way, that L-R is already high-pass filtered anyway, by the notch filter which reduces the carrier after AM modulation. So low frequencies are already imperfect anyway.

[gf]

Another interesting design is a keyed stereo decoder that uses time-multiplexing:
...
It is noted that the shorter the keying pulses, the greater the channel separation and the more accurately the A and B signals are reproduced. To achieve >30 dB separation and <1% distortion over 0.3–5 kHz, the control pulse width must not exceed 5 µs.

In other words, the composite signal is sampled at the positive and negative peaks of the subcarrier sine wave. Then, L and R are reconstructed with a brickwall lowpass filter. I tried that, and it does indeed seem to work very well.

However, I don't find it well suited for a DSP implementation because it requires a very high sample rate for precise timing of the sampling points if the sample rate and subcarrier frequency are unrelated. Perhaps interpolation could mitigate this issue, but since the subcarrier needs to be recovered anyway, downconversion via multiplying by the recovered complex subcarrier is easier.

[radiolistener]

Could you provide an example for such device?

It can be any SDR with transmit capability that supports an IQ baseband bandwidth of at least 384 kHz — for example, a PlutoSDR, HackRF One, or similar hardware.

As an alternative, this could also be tested with your receiver by using a recording of the stereo composite itself instead of the IQ FM signal. In that case, you would need any playback device capable of reproducing a 50 kHz bandwidth signal from a digital recording and feeding it into the appropriate test point of the receiver. However, such playback devices may be hard to find in practice, since typical audio interfaces include a low-pass filter with a cutoff around 20 kHz.

Here is actual spectrum of the test signal from this file.

[radiolistener]

Ideally, if the composite signal is exactly k1*(L+R)/2 + k2*subcarrier*(L-R)/2 + 0.19953 * subcarrier, where subcarrier has a peak amplitude of 1, then the DC offset (after down-conversion by multiplication with exp(-2i*pi*f_subcarrier*t)) should be a fixed constant value of 0.19953/2 = 0.099763, which you could subtract.

According to the book, the exact coefficients are k1=k2=0.8. However, it states that these factors are applied directly to the L+R and L−R signals, rather than to the DSB/AM component. The remaining portion of the amplitude is allocated to the subcarrier, which is given as 0.2 according to the book.

Thus, the expression takes the form 0.8*M + ((0.8*S) .* subcarrier) + 0.2*subcarrier
where
\$M = \frac{L + R}{2}\$;  \$S = \frac{L - R}{2}\$

With this proportion, the resulting composite has a peak amplitude of exactly 1, which I have verified — so in terms of scaling factors, everything checks out.

The only practical question is how the subcarrier is attenuated. One option is to generate the subcarrier separately, attenuate it precisely to 0.2, and then add it to the DSB-SC: (0.8*S) .* subcarrier + 0.2*subcarrier. This is spectrally clean and does not affect the signal shape.

Alternatively, as described in the book, transmitters can use AM modulation (DSB-C): (0.8*S + 1) .* subcarrier, in which case the subcarrier does not need to be added separately. However, this approach requires spectral shaping of the signal, which according to the book is implemented using an RC network with Q = 100 ± 5, essentially an IIR filter.

The first approach is more in line with a DSP implementation, while the second corresponds to the design philosophy of traditional analog transmitters. However, given that transmitters are designed to employ such a filter, it may be more appropriate to adopt the same approach in DSP as well, in order to ensure proper reconstruction of the low-frequency components of the signal.

However, in the case of imperfections, a high pass filter is likely the only option for removing an unknown DC component. You could use a very low cutoff, say 1 Hz. Note, by the way, that L-R is already high-pass filtered anyway, by the notch filter which reduces the carrier after AM modulation. So low frequencies are already imperfect anyway.

Yes, in general this approach does work. However, I have encountered a subtle issue: it produces a faint, periodically varying hum. As I mentioned earlier, I provided an example recording of the decoded signal where this hum is present and you can hear it in pauses. I am not entirely sure of its origin, but it appears to be linked to the composite signal’s sampling rate: at Fs = 192 kHz the hum is noticeable, whereas at Fs = 187.5 kHz it disappears. My current suspicion is that this may be related to noise introduced by the DSP PLL when operating at sampling rates that are not an integer multiple of 31.25 kHz, but I cannot rule out the possibility of a mistake somewhere in my implementation. It is also possible that this background hum is present in the generated composite itself. 

[gf]

The remaining portion of the amplitude is allocated to the subcarrier, which is given as 0.2 according to the book.

Obviously, the 0.2 stands for the -14dB reduced carrier (exact value would be 0.19953).

The only practical question is how the subcarrier is attenuated. One option is to generate the subcarrier separately, attenuate it precisely to 0.2, and then add it to the DSB-SC: (0.8*S) .* subcarrier + 0.2*subcarrier. This is spectrally clean and does not affect the signal shape.

Alternatively, as described in the book, transmitters can use AM modulation (DSB-C): (0.8*S + 1) .* subcarrier, in which case the subcarrier does not need to be added separately. However, this approach requires spectral shaping of the signal, which according to the book is implemented using an RC network with Q = 100 ± 5, essentially an IIR filter.

The first approach is more in line with a DSP implementation, while the second corresponds to the design philosophy of traditional analog transmitters. However, given that transmitters are designed to employ such a filter, it may be more appropriate to adopt the same approach in DSP as well, in order to ensure proper reconstruction of the low-frequency components of the signal.

I wonder if it is possibly easier for a PLL to lock if there is a gap around the subcarrier?

However, in the case of imperfections, a high pass filter is likely the only option for removing an unknown DC component. You could use a very low cutoff, say 1 Hz. Note, by the way, that L-R is already high-pass filtered anyway, by the notch filter which reduces the carrier after AM modulation. So low frequencies are already imperfect anyway.

Yes, in general this approach does work. However, I have encountered a subtle issue: it produces a faint, periodically varying hum. As I mentioned earlier, I provided an example recording of the decoded signal where this hum is present and you can hear it in pauses. I am not entirely sure of its origin, but it appears to be linked to the composite signal’s sampling rate: at Fs = 192 kHz the hum is noticeable, whereas at Fs = 187.5 kHz it disappears. My current suspicion is that this may be related to noise introduced by the DSP PLL when operating at sampling rates that are not an integer multiple of 31.25 kHz, but I cannot rule out the possibility of a mistake somewhere in my implementation. It is also possible that this background hum is present in the generated composite itself. 

Hmmm, there's not much the encoder/modulator can do wrong. 192 kSa/s should be fine. And the effect of high-pass filtering the L-R channel should only reduce the stereo width at low frequencies - it is not supposed to introduce (nonlinear) distortion.

To rule out the encoder/decoder, try demodulating without using the PLL-recovered subcarrier. Instead, use the same complex subcarrier, exp(-2i*pi*f_(subcarrier)*t), whose real part was used by the encoder for the modulation.

Btw, how do you demodulate/decode?
  L = composite .* (1 + 2 * recovered_complex_subcarrier)
  R = composite .* (1 - 2 * recovered_complex_subcarrier)
  % where recovered_complex_subcarrier is a complex sine wave exp(-2i*pi*f*t)
  % with correct frequency and phase
  Finally, DC-removal + de-emphasis + lowpass filterting of  L and R with
  a brickwall, cutting off steeply with a 15 kHz ... 16.25 kHz transition band)
?

I rather suspect the PLL. If it does not recover a clean carrier, the down-conversion can well produce intermodulation products.
Check the recovered subcarrier with a FFT (after the PLL has settled).

[gf]

To rule out the encoder/decoder, try demodulating without using the PLL-recovered subcarrier. Instead, use the same complex subcarrier, exp(-2i*pi*f_(subcarrier)*t), whose real part was used by the encoder for the modulation.

Btw, how do you demodulate/decode?
  L = composite .* (1 + 2 * recovered_complex_subcarrier)
  R = composite .* (1 - 2 * recovered_complex_subcarrier)
  % where recovered_complex_subcarrier is a complex sine wave exp(-2i*pi*f*t)
  % with correct frequency and phase
  Finally, DC-removal + de-emphasis + lowpass filterting of  L and R with
  a brickwall, cutting off steeply with a 15 kHz ... 16.25 kHz transition band)
?

The following works nicely for me:
[ pre/de-emphasis not yet implemented ]

Code: [Select]


pkg load signal

% read stereo audio file, the file must alreay
% have sufficient sample rate for composite
[x,fs] = audioread("/tmp/test1.wav");
L = x(:,1)';
R = x(:,2)';

% 15 kHz brickwall lowpass
[n,Wn,beta,ftype] = kaiserord([14500 15000],[1 0],[0.0001 0.0001],fs);
n += mod(n,2);
lp15k = fir1(n,Wn,ftype,kaiser(n+1,beta),"noscale");

% limit audio bandwith
L = fftfilt(lp15k,L);
R = fftfilt(lp15k,R);

% generate subcarrier
N = length(L)
fsubc = 31250
subc = exp(-2i*pi*fsubc*[0:length(L)-1]/fs);

% stereo encoder
composite = 0.8*(L+R)/2 + 0.8*real(subc).*(L-R)/2 + 0.2*real(subc);
audiowrite("/tmp/composite.wav", composite', fs);

% stereo decoder
LL = composite .* (1 + 2 * subc);
RR = composite .* (1 - 2 * subc);
LL -= mean(LL);
RR -= mean(RR);
LL = fftfilt(lp15k,LL);
RR = fftfilt(lp15k,RR);
audiowrite("/tmp/decoded.wav", real([LL;RR])', fs);


Test signal is the first one from here: https://pixabay.com/music/search/30%20sec/

[gf]

Generated composite signal.

[gf]

Stereo decoded from the generated composite, using the original subcarrier from the modulator.

[gf]

I also tried playing with a PLL, based on this document [ I did some modifications, though ]:
https://s3.amazonaws.com/embeddedrelated/user/113580/digital%20plls%20part%203_68269.pdf
My sample rate (176400) is not related to the subcarrier frequency.

[ Btw, in Octave, the PLL code runs horribly slow, as it is necessary to loop over the samples instead of working with matrices. ]

Attached is the L/R decoded from my generated composite, using a PLL-recoverd subcarrier for the demodulation.
With an intentionally wrong initial frequency (offset of 100 ppm), it locks within a few seconds. It doesn't seem to work that badly.

[ However, so far I had no luck when I tried to lock it to Njk's composite signal and to decode it. The result is significantly distorted. ]

[radiolistener]

[ Btw, in Octave, the code runs horribly slow, as it is necessary to loop over the samples instead of working with matrices. ]

For such real-time-like processing, a common approach is to implement the heavy part in C++ as a MEX function, which can then be called directly from Octave almost as fast as native code. The only con is that each time you modify CPP file, you need to recompile it - just execute command like:

Code: [Select]

mkoctfile --mex PLL_MEX.cpp

it will generate PLL_MEX.mex and you can call it from octave as PLL_MEX(...).

For example, here is one of my old attempts to implement PLL (I'm not very familiar with PLL math, so it may be mistaken). Currently I'm using another one based on liquid-dsp library, but its less readable.

1)The MEX extension code

Code: [Select]

// linux: mkoctfile --mex PLL_MEX.cpp
#include "mex.h"
#include <cmath>

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
    // Inputs:
    // prhs[0]  = samples (double vector)
    // prhs[1]  = w0
    // prhs[2]  = min_w
    // prhs[3]  = max_w
    // prhs[4]  = alpha
    // prhs[5]  = beta
    // prhs[6]  = phaseAdj
    // prhs[7]  = lockAlpha
    // prhs[8]  = lockOneMinusAlpha
    // prhs[9]  = lockThreshold
    // prhs[10] = phase
    // prhs[11] = phaseErrorAvg

    double *samples = mxGetPr(prhs[0]);
    size_t N = mxGetNumberOfElements(prhs[0]);

    double w0 = mxGetScalar(prhs[1]);
    double min_w = mxGetScalar(prhs[2]);
    double max_w = mxGetScalar(prhs[3]);
    double alpha = mxGetScalar(prhs[4]);
    double beta  = mxGetScalar(prhs[5]);
    double phaseAdj = mxGetScalar(prhs[6]);
    double lockAlpha = mxGetScalar(prhs[7]);
    double lockOneMinusAlpha = mxGetScalar(prhs[8]);
    double lockThreshold = mxGetScalar(prhs[9]);

    double phase = mxGetScalar(prhs[10]);
    double phaseErrorAvg = mxGetScalar(prhs[11]);

    // Outputs
    plhs[0] = mxCreateDoubleMatrix(N,1,mxCOMPLEX); // vco
    plhs[1] = mxCreateDoubleMatrix(N,1,mxREAL);    // isLocked
    plhs[2] = mxCreateDoubleMatrix(N,1,mxREAL);    // adjustedPhase
    plhs[3] = mxCreateDoubleMatrix(N,1,mxREAL);    // phaseErrorAvg
    plhs[4] = mxCreateDoubleMatrix(1,3,mxREAL);    // updated state: phase, w0, phaseErrorAvg

    double *vco_re = mxGetPr(plhs[0]);
    double *vco_im = mxGetPi(plhs[0]);
    double *locked = mxGetPr(plhs[1]);
    double *adjPhase = mxGetPr(plhs[2]);
    double *peAvg = mxGetPr(plhs[3]);
    double *stateOut = mxGetPr(plhs[4]);

    for (size_t n=0; n<N; n++) {
        // generate VCO
        double c = cos(phase);
        double s = sin(phase);
        double x = samples[n];

        // Costas phase detector
        // I = input * cos(VCO)
        // Q = input * (-sin(VCO))
        double I = x * c;
        double Q = x * (-s);

        // phase error = I*Q (Costas detector)
        double phaseError = I * Q;

        // loop filter
        w0 += beta * phaseError;
        if (w0 > max_w) w0 = max_w;
        if (w0 < min_w) w0 = min_w;

        // exponential moving average of squared phase error
        phaseErrorAvg = lockOneMinusAlpha * phaseErrorAvg + lockAlpha * phaseError * phaseError;

        // update phase
        phase += w0 + alpha * phaseError;
        phase = fmod(phase, 2*M_PI);
        double adjP = phase + phaseAdj;

        // lock detection
        int isLock = (phaseErrorAvg < lockThreshold) ? 1 : 0;

        // store outputs
        vco_re[n] = c;
        vco_im[n] = s;
        adjPhase[n] = adjP;
        peAvg[n] = phaseErrorAvg;
        locked[n] = isLock;
    }

    // Return updated state
    stateOut[0] = phase;
    stateOut[1] = w0;
    stateOut[2] = phaseErrorAvg;
}

2) Usage:

Code: [Select]

                Fs = ...;
                f0 = 31250;  % target frequency Hz
                df0=20;       % frequency range Hz
                bw=10;       % bandwidth Hz
                zeta=0.707;
                lockTime=0.5; % seconds
                lockThreshold=1.0;

            w0 = 2*pi*f0 / Fs;
            min_w = (2*pi*f0 - df0) / Fs;
            max_w = (2*pi*f0 + df0) / Fs;
            alpha = 2*zeta * 2*pi*bw / Fs;
            beta  = alpha^2 / (4.0 * zeta^2);
            phaseAdj = -1.75;
            lockAlpha = 1.0 - exp(-1.0 / (Fs * lockTime));
            lockOneMinusAlpha = 1.0 - lockAlpha;
           
            phase = 0;
            adjustedPhase = 0;
            phaseErrorAvg = 0;
            isLocked = 0;


                % call MEX
                [vco, isLocked, adjustedPhase, phaseErrorAvg, state] = PLL_MEX(...
                    samples, w0, min_w, max_w, alpha, beta, ...
                    phaseAdj, lockAlpha, lockOneMinusAlpha, ...
                    lockThreshold, phase, phaseErrorAvg);


[radiolistener]

  L = composite .* (1 + 2 * recovered_complex_subcarrier)
  R = composite .* (1 - 2 * recovered_complex_subcarrier)

Cool, nice approach, I did it separately:

Code: [Select]

mixed = composite .* _osc3125;
mixed *= 2;           % correct amplitude

% remove DC
fc_dc = 10; % 10-20 Hz cutoff for DC estimation
[b_dc, a_dc] = butter(1, fc_dc/(Fs/2), 'low'); 
dc_est = filter(b_dc, a_dc, mixed);   % low-frequency component ~DC
mixed -= dc_est;

% extract L-R
_sub_LR = filter(lpf_15k, 1, mixed);
_sub_LR /= 0.4;

% extract L+R
_add_LR = filter(lpf_15k, 1, composite);
_add_LR /= 0.4;

_L = (_add_LR + _sub_LR) / 2;
_R = (_add_LR - _sub_LR) / 2;


[tszaboo]

I was not referring to the absolute subcarrier amplitude. For proper stereo decoding it’s essential that the L+R and L–R paths maintain equal gain. The L+R component is straightforward, but the L–R signal first needs to be separated from the carrier and demodulated from DSB. This can be done with a synchronous detector using the PLL-locked carrier. However, because the carrier is present within DSB signal bandwidth, the result includes a DC offset and require HPF to remove it.

I’d prefer not to use a high-pass filter just to remove that DC. My idea was to try subtracting the carrier directly, by taking the PLL output and attenuating it by 14 dB, then substract it from composite. The issue is that the actual level of the carrier may vary between transmitters, so a fixed 14 dB attenuation would likely leave a residual offset due to this mismatch.

Any idea on how to extract L-R from DSB which include carrier but with no use of HPF?

What happens mathematically, if you demodulate the signal 31.25 kHz/2 offset from the carrier frequency? It's too late right now to think clearly, as it's almost midnight, but I think using the mirror image (negative frequency) of the L-R signal and mixing it with the L+R might be somehow simple.

[Njk]

Looks like it's implemented in the gqrx somehow (starting from v.2.4)
https://github.com/gqrx-sdr/gqrx/blob/master/src/dsp/stereo_demod.cpp

[gf]

Looks like they do bandpass-filtering, before feeding into the PLL (which I did not do in my experiments - I directly tried to lock to the composite), This certainly improves locking and far-from-carrier noise (it won't improve close-in phase noise, though). Maybe I should try this as well.

[radiolistener]

I tried installing gqrx. The interface is a bit unusual, but I managed to feed it a file. It runs rather slowly—unfortunately my machine doesn’t have enough performance for smooth playback without dropouts—but it does decode OIRT stereo. However, it seems to swap the L and R channels, and the stereo detection doesn’t always trigger. It’s unclear whether this is due to immature code or a misconfiguration on my part. But at a glance it appears to achieve good channel separation, so it’s worth studying the code more carefully.

[gf]

Can the gqrx stereo decoder demodulate Njk's recorderd composite w/o distortion?

[radiolistener]

Can the gqrx stereo decoder demodulate Njk's recorderd composite w/o distortion?

I removed it because with all its libraries and dependencies it was taking up quite a lot of disk space. However, I’m planning to extract just the OIRT/CCIR stereo decoding part into a separate test program for experimentation and will try it.