I experimented with several different phase detectors. Some of them can work reasonably well for CCIR stereo, but OIRT turned out to be much more sensitive to accurate carrier phase recovery. Even relatively small phase errors noticeably degrade the L-R demodulation.
Because of that, I currently use a conjugate mixer to shift the pilot to DC and then calculate the phase error with atan2(). In my experiments this produced the best results for OIRT stability and stereo separation.
You may be right that I am missing some simpler or more lightweight approach, but so far all alternatives I tested gave worse results on real OIRT recordings, especially regarding phase accuracy.
Here is my current PLL implementation in C++ MEX for Octave.
In its current form, it only shows the PLL acquisition behavior for a clean synthetic sine wave. However, you can replace it with loading an audio file containing a real carrier signal for testing — examples of such file loading are already present in the code as commented sections.
pll_process_mex.cpp:
// linux: mkoctfile --mex pll_process_mex.cpp
#include "mex.h"
#include <cmath>
#include <complex>
typedef struct {
float phase; // radians
double freq; // radians/step
float err;
std::complex<float> lo;
float lock_avg;
bool is_lock;
float alpha;
float beta;
double min_freq;
double max_freq;
} pll_t;
static inline float mod_2pi(float x) {
if (x > M_PI) return x - (2.0 * M_PI);
else if (x < -M_PI) return x + (2.0 * M_PI);
else return x;
}
static void pll_init(pll_t &pll, double loop_bw, double min_freq, double max_freq) {
if (loop_bw < 0) {
throw std::out_of_range("pll_init: invalid bandwidth. Must be >= 0.");
}
pll.phase = 0.0f;
pll.freq = 0.0f;
pll.err = 0.0f;
pll.lo = std::exp(std::complex<float>(0.0f, pll.phase));
pll.lock_avg = 0.0;
pll.is_lock = false;
pll.max_freq = max_freq;
pll.min_freq = min_freq;
// update gains
double damping = sqrt(2.0) / 2.0;
double nfactor = (1.0 + 2.0 * damping * loop_bw + loop_bw * loop_bw);
pll.alpha = (4.0 * damping * loop_bw) / nfactor;
pll.beta = (4.0 * loop_bw * loop_bw) / nfactor;
}
static inline void pll_update(pll_t &pll, std::complex<float> input) {
// oscillator
float sin, cos;
sincosf(pll.phase, &sin, &cos);
pll.lo = std::complex<float>(cos, sin);
// phase detector
std::complex<float> pdiff = input * std::conj(pll.lo);
pll.err = atan2f(pdiff.imag(), pdiff.real());
pll.freq += pll.beta * pll.err; // freq update
pll.freq = std::min(std::max(pll.freq, pll.min_freq), pll.max_freq); // freq clamp
pll.phase += pll.freq + pll.alpha * pll.err;// phase update
pll.phase = mod_2pi(pll.phase); // phase wrap
// lock detect
float detect = std::abs(pll.err); // abs(err)
detect = 1.0f - std::min(detect, 1.0f); // clamp 0..1 and invert
pll.lock_avg += pll.alpha * (detect - pll.lock_avg); // IIR LPF
if (pll.is_lock && pll.lock_avg < 0.4f) {
pll.is_lock = false;
fprintf(stderr, "PLL loss\n");
} else if (!pll.is_lock && pll.lock_avg > 0.6f) {
pll.is_lock = true;
fprintf(stderr, "PLL lock\n");
}
}
//----------------------------------------------------------------------
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
if (nrhs != 4) {
mexErrMsgTxt("Usage: [output, err, lock] = pll_process_mex(input, loop_bw, min_freq, max_freq)");
}
const mxArray *in = prhs[0];
double *x_re = mxGetPr(in);
double *x_im = mxGetPi(in);
if (!x_im) {
mexErrMsgTxt("input must be complex");
}
double loop_bw = mxGetScalar(prhs[1]);
double min_freq = mxGetScalar(prhs[2]);
double max_freq = mxGetScalar(prhs[3]);
mwSize N = mxGetNumberOfElements(in);
// Create output arrays
plhs[0] = mxCreateDoubleMatrix(N, 1, mxCOMPLEX);
double *y_re = mxGetPr(plhs[0]); // output
double *y_im = mxGetPi(plhs[0]);
plhs[1] = mxCreateDoubleMatrix(N, 1, mxREAL); // err
double *err_re = mxGetPr(plhs[1]);
plhs[2] = mxCreateDoubleMatrix(N, 1, mxREAL); // lock
double *lock_re = mxGetPr(plhs[2]);
// init
pll_t pll;
pll_init(pll, loop_bw, min_freq, max_freq);
for (mwSize n = 0; n < N; n++) {
pll_update(pll, std::complex<float>(x_re[n], x_im[n]));
y_re[n] = pll.lo.real();
y_im[n] = pll.lo.imag();
err_re[n] = pll.err;
//lock_re[n]= pll.is_lock;
lock_re[n]= pll.lock_avg;
}
}
And here is Octave test code for PLL:
pkg load signal;
fs = 384000; % sample rate
f0 = 31250; % OIRT pilot
N = fs; % 1 sec test
t = (0:N-1)'/fs;
% --- test signal: pure pilot tone ---
phi0 = pi/3; % initial phase of pilot in radians
pilot = cos( (2*pi*f0*t + phi0));
%[pilot, fs] = audioread('/home/pi/tmp-oirt/50k/test-oirt-1-lpf96k.wav', [1 N]); % read file
%[pilot, fs] = audioread('../../DATA/test-oirt-composite-50us.flac', [1 N]); % read file
if size(pilot, 2) ~= 1
error('Input file must be mono composite');
end
N = length(pilot);
%---UTILS---------------------------------------------------------------
% Kaiser beta estimator
function bta = firdes_kaiser_beta(A)
% A - sidelobe suppression in dB
if (A >= 50)
bta = 0.1102 * (A - 8.7);
elseif (A < 50 && A > 21)
bta = 0.5842 * (A - 21)^0.4 + 0.07886 * (A - 21);
else
bta = 0;
end
end
% FIR Kaiser ntaps estimator
function ntaps = firdes_kaiser_ntaps(fs, df, A)
dw = 2*pi * df / fs;
if dw <= 0
error('df must be > 0 Hz');
end
ntaps = max(1, ceil((A-7.95)/(2.285*dw)));
end
% FIR real LPF (FIR Type I)
function taps = firdes_real_lpf(gain, fs, fc, tw, A)
N = firdes_kaiser_ntaps(fs, tw, A);
% force odd number of taps for Type I FIR
% (odd length, symmetric impulse response, LPF/BPF)
if mod(N,2)==0
N = N + 1;
end
M = floor((N-1)/2);
% ideal sinc LPF impulse response
wc = 2*pi*fc/fs;
taps = zeros(1, N);
n=-M:M;
idx = (n ~= 0);
taps(idx) = sin(wc*n(idx)) ./ (pi*n(idx)); % h(n) = sin(wc*n)/(pi*n)
taps(~idx) = wc/pi; % h(0) = wc/pi
% apply Kaiser window
taps = taps .* kaiser(N, firdes_kaiser_beta(A)).';
% normalize DC gain
fmax = taps(M+1) + 2*sum(taps(M+2:end));
gain = gain / fmax;
taps = taps * gain;
end
% FIR complex BPF
function taps = firdes_complex_bpf(gain, fs, f1, f2, tw, A)
% prototype lpf
lpf_taps = firdes_real_lpf(gain, fs, (f2-f1)/2, tw, A);
fc = (f1 + f2)/2; % center freq [Hz]
wc = 2*pi*fc/fs; % center freq [rad/sample]
N = length(lpf_taps);
if mod(N,2)==0
error('firdes_complex_bpf: odd taps LPF expected');
end
M = floor((N-1)/2);
n = -M:M;
taps = lpf_taps .* exp(1j * wc * n); % shift to wc
end
% plot signal FFT with Kaiser window
% A - sidelobe suppression in dB
function plot_fft_kaiser(x, fs, A, text)
N = length(x);
% truncate to NPOT
NPOT = 2^floor(log2(N));
NPOT = min(NPOT, 2^22);
if N ~= NPOT
N = NPOT;
x = x(1:N);
end
fprintf('FFT: %g kpts\n', N/1024);
xw = x(:) .* kaiser(N, firdes_kaiser_beta(A)); % apply window
x_fft = fftshift(fft(xw));
mag = abs(x_fft);
mag = mag / max(mag + eps); % normalize
mag_db = 20.0*log10(mag + eps);
f = ((0:N-1) - floor(N/2)) * fs/N;
% autoscale freq axis
limits = [1e7 1e4 0]; scales = [1e6 1e3 1]; snames = {'M','k',''};
i = find(max(abs(f)) >= limits, 1, 'first');
if isempty(i), i = numel(scales); end
f = f / scales(i);
funit = sprintf('%sHz', snames{i});
% plot
figure;
plot(f, mag_db);
grid on; grid minor;
set(gca, 'gridalpha', 0.2); set(gca, 'minorgridalpha', 0.05);
xlabel(funit);
ylabel('dB');
xlim([min(f) max(f)]);
ylim([min(mag_db)-6 max(mag_db)+6]);
title(sprintf('%s\nFFT %g kpts, fs=%g kS/s', text, N/1024, fs/1000));
end
% Ensure MEX is built and up-to-date
function ensure_mex_uptodate(src_name)
[pathstr, name, ext] = fileparts(src_name);
if isempty(ext)
error('Source file must include extension');
end
mex_name = [name, '.', mexext()]; %mex_name = fullfile(pathstr, [name, '.', mexext()]);
src = dir(src_name);
mex = dir(mex_name);
if isempty(src)
error('Source file not found: %s', src_name);
end
if isempty(mex) || src.datenum > mex.datenum % Rebuild if needed
fprintf('Compile %s...\n', src_name);
status = system(sprintf('mkoctfile --mex %s', src_name));
if status ~= 0
error('MEX compile failed');
end
end
end
%-----------------------------------------------------------------------
% --- initialize PLL/BPF ---
pll_bw = 5;
pll_minf = f0-50;
pll_maxf = f0+50;
bpf = firdes_complex_bpf(1.0, fs, pll_minf, pll_maxf, 5000.0, 60);
bpf_delay = (length(bpf) - 1) / 2; % bpf_delay = mean(grpdelay(bpf));
fprintf('FIR BPF: %d taps, delay = %f samples\n', length(bpf), bpf_delay);
bpf_delay = floor(bpf_delay);
%freqz(bpf, 1, fs, 2^16);
% apply BPF
pilotc = complex(pilot*2.0, 0); % to complex with amp compensation (mirror will be removed with BPF)
pilotc = filter(bpf, 1, pilotc);
fprintf('bpf done\n');
% fix for BPF delay
pilotc = [pilotc(1+bpf_delay:end); zeros(bpf_delay, 1)];
ensure_mex_uptodate('pll_process_mex.cpp');
[pll_lo, err, lock] = pll_process_mex(pilotc, pll_bw * 2.0*pi/fs, pll_minf * 2.0*pi/fs, pll_maxf * 2.0*pi/fs);
fprintf('pll done\n');
% plot FFT
lock_idx = find(lock > 0.6, 1);
if ~isempty(lock_idx)
plot_fft_kaiser(pll_lo(lock_idx:end), fs, 100, ...
sprintf('loop\\_bw=%g, minf=%g, maxf=%g', pll_bw, pll_minf, pll_maxf));
end
aerr = real(pilotc) - pilot;
% --- plot Waveform ---
%t = 1:length(pilotc); % sample number
t = (1:length(pilotc)) / fs; % second
figure;
hold on;
%plot(t, real(pilotc), 'DisplayName', 'pilotc');
%plot(t, real(pll_lo), 'DisplayName', 'pll\_lo');
plot(t, aerr, 'DisplayName', 'amplitude-err');
plot(t, err, 'DisplayName', 'phase-err');
plot(t, lock, 'DisplayName', 'lock');
grid on; grid minor;
set (gca, 'gridalpha', 0.2); set(gca, 'minorgridalpha', 0.05);
xlim([min(t) max(t)]);
%xlim([95830 95863]); ylim([-1.1 1.1]);
ylim([-0.6 2.6]);
xlabel('sec');
legend();
title(sprintf('loop\\_bw=%g, minf=%g, maxf=%g', pll_bw, pll_minf, pll_maxf));
Home
Help
Search
Login
Register
