
composite = 0.8*add_LR + (0.8*sub_LR .* osc3125) + 0.2*osc3125;
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).
sox test-oirt-wfm-50us.flac test-oirt-wfm-50us.wav
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)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?
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?
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.
Could you provide an example for such device?
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*fsubcarrier*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.
The remaining portion of the amplitude is allocated to the subcarrier, which is given as 0.2 according to the book.
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.
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)
?
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);
[ Btw, in Octave, the code runs horribly slow, as it is necessary to loop over the samples instead of working with matrices. ]
mkoctfile --mex PLL_MEX.cpp
// 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;
}
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);
L = composite .* (1 + 2 * recovered_complex_subcarrier)
R = composite .* (1 - 2 * recovered_complex_subcarrier)
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;
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?
Can the gqrx stereo decoder demodulate Njk's recorderd composite w/o distortion?