EEVblog® Electronics Community Forum
Electronics => RF, Microwave, Ham Radio => Topic started by: radiolistener on August 24, 2025, 02:51:28 pm
-
As you probably know, there are two standards for WFM stereo broadcasting - one applied in the OIRT band (65.9 - 74 MHz) and the other in the CCIR band (87.5 - 108.0 MHz).
OIRT uses 31.25 kHz subcarrier to transmit L-R component, with so called "polar modulation" and partially suppressed subcarrier.
CCIR uses 38 kHz subcarrier to transmit L-R component, with DSB-SC and pilot-tone at 19 kHz in order to recover subcarrier.
I am working on adding OIRT FM stereo decoder support to my receiver. Unfortunately, there is no longer any stereo broadcasting in the OIRT band in my region - only mono transmission remains. I also don’t have access to IQ recordings of actual OIRT stereo broadcasts for testing, so I have to rely on old magazine articles and schematics of stereo decoders.
One point I find confusing is the terminology: in the OIRT stereo composite signal, the modulation is often described as "polar modulation". It is not entirely clear to me what exactly is meant by this.
From what I understand, the stereo composite modulation in both CCIR and OIRT systems is essentially DSB. The main difference is that in CCIR, the 38 kHz subcarrier is completely suppressed and recovered on the receiver side by doubling the transmitted 19 kHz pilot tone, whereas in OIRT stereo-composite there is no pilot tone - instead, the 31.25 kHz subcarrier is only partially suppressed (about -14 dB), allowing it to be recovered directly on the receiver side.
Apart from this subcarrier recovery method, it seems to me that the modulation principle is the same in both systems, so the same DSB demodulator could be used for both standards, with the only difference being the subcarrier recovery mechanism.
Am I understanding this correctly?
If anyone happens to have a recording of an FM broadcast in OIRT stereo mode (65.9 - 74 MHz band), I would greatly appreciate it. Ideally this could be provided as a WAV file containing an IQ IF stream, or alternatively as a direct recording of the composite signal from the FM detector output.
-
Here's what Google Gemini told me:
https://g.co/gemini/share/2583583ca5ad
I can't assess correctness, but it sounds plausible.
-
Here's what Google Gemini told me:
https://g.co/gemini/share/2583583ca5ad
I can't assess correctness, but it sounds plausible.
It seems the Gemini AI response introduces some confusion here. It states:
In the OIRT system, polar modulation is used to modulate the main carrier, not the subcarrier.
which is clearly incorrect. In fact, polar modulation in the OIRT stereo system refers to the modulation of the stereo subcarrier in the composite signal, not the main FM carrier itself.
If I understand correctly, the mathematical difference can be described as follows (Octave code), we have:
Fs = 384000; % sample rate
N = length(L) % sample count
t = (0:N-1).' / Fs; % time vector
% Build sum and difference
add_LR = (L + R) / 2; % sum channel
sub_LR = (L - R) / 2; % difference channel
Then:
- CCIR uses this:
pilot = 19000; % pilot tone frequency 19 kHz
subcf = pilot*2; % sub-carrier frequency 38 kHz
pilot = cos(2*pi*pilotf*t); % pilot tone 19 kHz
c38 = -sin(2 * 2*pi*pilotf*t); % carrier 2*19=38 kHz
dsbsc = sub_LR .* c38;
% Composite baseband signal (L+R + pilot + DSB-SC)
composite = 0.45*add_LR + 0.1*pilot + 0.45*dsbsc;
While OIRT uses something like this:
subcf = 31250; % sub-carrier frequency 31.25 kHz
c3125 = sign(cos(2*pi*subcf*t)); % squarewave subcarrier
polar = 0.2*c3125 + sub_LR .* c3125;
% Apply bandpass filter with center frequency 31.25 kHz ±15 kHz
polar = filtfilt(bpf_31250, 1, polar);
% Composite baseband signal (L+R + polar)
composite = 0.45*add_LR + 0.45*polar;
And since it uses squarewave carrier with following BPF, it is actually equals to:
c3125 = cos(2*pi*subcf*t); % sine wave subcarrier
polar = 0.2*c3125 + sub_LR .* c3125;
% Composite baseband signal (L+R + polar)
composite = 0.45*add_LR + 0.45*polar;
which is more easy and clean, because don't involve BPF to remove harmonics added with square wave carrier.
As result, in the OIRT case, the square wave subcarrier generates harmonics, which must then be removed by a band-pass filter to keep the spectrum within ±15 kHz. After filtering, this is essentially equivalent to using a DSB-SC scheme with a partially suppressed carrier. By contrast, the CCIR system is simply DSB-SC with the carrier fully suppressed and reconstructed from the 19 kHz pilot tone.
So in effect, the only difference is that OIRT leaves a residual carrier (at about -14 dB) while CCIR completely suppresses it. Also it means that you can use polar demodulator to decode both - OIRT and CCIR system. The only difference is a subcarrier recovery scheme.
Am I understanding this correctly?
I do not have access to the real OIRT stereo signal recording to verify my understanding of the stereo composite modulation in this system, so I have to rely on fragmented descriptions and circuit diagrams from old radio magazine publications from 1970...
Here is one of the schematics published in Radio magazine, issue No. 3, 1974. According to the circuit, transistor T1 amplifies the signal, T2 recovers the 31.25 kHz subcarrier, then the low-frequency part (L+R) is extracted, while the high-frequency part containing (L–R) is amplified by T3 (L3 is tuned at 31.25 kHz, Q=4.9, R16 is used to decrease Q-factor) and demodulated using a diode bridge. Finally, the L and R channels are formed as:
L = (add_LR + sub_LR) / 2;
R = (add_LR - sub_LR) / 2;
-
Here's what Google Gemini told me:
https://g.co/gemini/share/2583583ca5ad
I can't assess correctness, but it sounds plausible.
It seems the Gemini AI response introduces some confusion here. It states:
In the OIRT system, polar modulation is used to modulate the main carrier, not the subcarrier.
which is clearly incorrect. In fact, polar modulation in the OIRT stereo system refers to the modulation of the stereo subcarrier in the composite signal, not the main FM carrier itself.
I really don't known. However, polar modulation makes only sense to me in a context where a single carrier is modulated with two modulating signals. And this applies rather to the main carrier. Since the sub-carrier is only modulated with a single modulating signal (L-R), I see no need for polar modulation for the sub-carrier at all. What would be the 2nd modulating signal?
Of course, instead of polar modulation, L+R and the modulated sub-carrier could also be added, and then FM-modulated onto the main carrier (like CCIR). But then no polar modulation would be in the play at all. :-//
-
Perhaps someone could help by performing a practical test: the idea is to transmit a generated composite signal with a low-power transmitter and check whether OIRT stereo receivers produce correct stereo decoding. The test signal could be prepared either as a baseband composite signal for direct input into an FM modulator, or as a fully FM-modulated IQ recording centered at 0 Hz, which would then only need to be shifted to the desired frequency.
I checked through my old radio receivers, but unfortunately all of them are equipped only with CCIR stereo decoders, so at the moment I have nothing suitable to test with.
If anyone has an old radio receiver with an OIRT stereo decoder, or access to a local OIRT stereo broadcast, I would greatly appreciate it if you could let me know.
-
After challenging Gemini and pointing out contradictions, it changed its mind several times and its most recent conclusion was
Yes, based on all technical sources, pure FM for the main carrier is the correct and only method used in both the OIRT and CCIR standards.
The "polar modulation" that has caused so much confusion does not refer to the modulation of the main carrier. It is a historical term used to describe a specific way that the stereo information is encoded and combined with the mono signal before it ever reaches the main transmitter.
Main Carrier vs. Composite Signal
Think of the process in two steps:
The Composite Baseband Signal: This is the complete audio signal, which contains all the necessary components for both mono and stereo reception.
It includes the L+R signal (the mono sound).
It includes the stereo subcarrier (the L-R signal modulated onto a subcarrier, either at 31.25 kHz for OIRT or 38 kHz for CCIR).
It might also include other data signals like RDS.
The Main FM Carrier: This is the powerful radio wave that is actually broadcast from the tower. The entire composite baseband signal from the first step is used to frequency-modulate (FM) this main carrier. The main carrier's frequency is shifted up and down according to the instantaneous amplitude of the composite signal.
In both the OIRT and CCIR systems, the main carrier's modulation is strictly FM. The difference between the two systems is in what makes up the composite signal, not in how that composite signal modulates the main carrier.
Based on my sources, the OIRT FM stereo subcarrier uses Amplitude Modulation (AM). Specifically, it uses Double Sideband Reduced Carrier (DSB-RC) AM.
The alternative, and incorrect, claim that it uses FM is a major point of confusion. Here is the reason why the AM explanation is technically sound and the FM explanation is not:
Synchronization: The primary difference between the OIRT and CCIR standards is how the receiver decodes the stereo signal. The OIRT system's 31.25 kHz subcarrier is intentionally transmitted at a low, "reduced" amplitude (about 14 dB below the main signal). A receiver's phase-locked loop (PLL) uses this partially transmitted carrier to synchronize with and demodulate the stereo signal.
Technical Impossibility of FM Suppression: The concept of "partially suppressing" a carrier only applies to amplitude-modulated signals. The amplitude of an FM signal is constant; all the information is encoded in its frequency changes. You cannot reduce its amplitude for synchronization without losing the information it's carrying.
The term "polar modulation" is a historical misnomer that has contributed to this confusion. In this context, it describes the OIRT system's method of combining stereo signals in a way that is different from the CCIR system, but it does not accurately describe the subcarrier's actual modulation. The subcarrier itself is still amplitude-modulated.
I don't have direct access to the original, official standards for the OIRT or CCIR systems. These are historical documents not available in my training data, which consists of publicly accessible online text. My information is based on secondary sources, such as technical summaries, forum discussions, and educational texts that describe the standards. This reliance on secondary sources is why my previous explanations contained inconsistencies.
So I think your assumption is correct. It's likely similar to CCIR, but w/o pilot tone, and with AM/DSB-RC instead of AM/DSB-SC for the subcarrier, and the term "polar modulation" is just used to create confusion.
Today, the term 'polar modulation' does not refer to a specific type of modulation, but rather to a particular power amplifier (PA) architecture in which a linear PA is emulated using a highly efficient non-linear PA. The same RF signal that can be generated using I/Q modulation and a linear PA can also be generated using a polar modulation PA, with better efficiency.
-
I came across this discussion on generating the OIRT stereo composite:
https://forums.stereotool.com/viewtopic.php?t=6233 (https://forums.stereotool.com/viewtopic.php?t=6233)
Unfortunately, there are few technical details and no example signals provided, but it seems that the authors of StereoTool have implemented a working version of polar stereo encoding for OIRT. If anyone here uses this tool, could you possibly create a small test wav file with a stereo composite signal in polar stereo format?
The discussion also references a standard that provides some details on deviation and the relative levels of the components in the composite stereo signals for both OIRT and CCIR systems:
https://www.itu.int/dms_pubrec/itu-r/rec/bs/R-REC-BS.450-3-200111-S!!PDF-E.pdf (https://www.itu.int/dms_pubrec/itu-r/rec/bs/R-REC-BS.450-3-200111-S!!PDF-E.pdf)
According to this document, the frequency deviation for OIRT stereo broadcast is 50 kHz instead of 75 kHz which is used in CCIR.
Also it mentioned OIRT stereo-composite proportions:
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).
2.1.2.3 The stereophonic multiplex signal is the sum of:
– the pre-emphasized signal, M;
– the sideband spectral components which are the product of amplitude-modulated unsuppressed carrier by a pre-emphasized signal S additionally transformed from the law \$\overline{K}(f)\$
– the sub-carrier with the amplitude reduced by 14 dB.
2.1.2.4 The amplitudes of the various components of the stereophonic multiplex signal, referred to the maximum amplitude of that signal (which corresponds to the maximum frequency deviation) are:
– signal M : maximum value 80% (A and B being equal, and in phase);
– signal S : maximum value 80% (A and B being equal but of opposite phase);
– reduced sub-carrier at 31.25 kHz; maximum residual amplitude 20%.
2.1.2.5 The frequency modulation is arranged in such a way that positive values of the multiplex signal correspond to a positive frequency deviation of the main carrier and negative values to negative frequency deviation.
A = left channel
B = right channel
M = (A+B) / 2
S = (A-B) / 2
Both M and S are pre-emphasised with pre-emphasis 50 us (standard pre-emphasis for Europe).
AM modulated signal (which consists of S message) is spectrum shaped with mentioned K(f) law in order to suppress carrier for -14 dB.
From the text it is not entirely clear whether the spectrum of the AM-modulated signal should be corrected according to K(f), or if the more correct approach is simply to combine the full DSB-SC amplitude with the carrier attenuated by -14 dB.
-
From the text it is not entirely clear whether the spectrum of the AM-modulated signal should be corrected according to K(f), or if the more correct approach is simply to combine the full DSB-SC amplitude with the carrier attenuated by -14 dB.
Yes it isn't clear. Interestingly, the transfer function K(f) is a highpass with ~0.8Hz cutoff and 14dB stopband attenuation. If the subcarrier amplitue is suposed to be represented by the DC component of the baseband modulating signal, then this filter - applied to the baseband signal - would provide the -14dB subcarrier attenuation.
-
If anyone has an old radio receiver with an OIRT stereo decoder, or access to a local OIRT stereo broadcast, I would greatly appreciate it if you could let me know.
Well, we still have 3 working stations of the OIRT band here in St. Petersburg. For one of them, 69.47 MHz, my old receiver indicates a stereo signal (just checked).
https://radiomap.eu/ru/sankt-peterburg
Perhaps you can use a WebSDR connection to view the signal
-
For one of them, 69.47 MHz, my old receiver indicates a stereo signal (just checked).
Perhaps you can use a WebSDR connection to view the signal
Unfortunately websdr doesn't support wideband FM, its limit - max 20 kHz, while OIRT stereo FM requires at least 96.25 kHz bandwidth. And there is no websdr receiver covering 69.47 MHz at that location, only ham bands on HF and VHF...
Could you please record the file with baseband IQ of this station with RTLSDR? You can do it with HDSDR (https://www.hdsdr.de/).
It is important to setup LO outside of station bandwidth - use LO=69.3 MHz or LO =69.7 MHz for recording, it helps to avoid DC spur which can affect signal quality. It's better to setup RTLSDR sample rate 960 kHz instead of 2.4 MHz to reduce file size, since I need just +-150 kHz around station frequency. Also please disable AGC and setup good RF gain manually, to avoid unwanted gain variations.
30-60 seconds record will be enough. If possible, please record it when station broadcasting some music or sound which allows to hear stereo effects and distinguish it from mono. If such music composition will allow to check proper left/right channel decoding it will be great.
You can share file on https://transfiles.ru/ (https://transfiles.ru/)
-
The problem is that it's old Soviet receiver, Leningrad RP-015 to be precise. It has no IQ outputs. The model is well known so you can find the schematic diagrams on the internet. There are no patented chips, it's made entirely of discrete components. There is a test point right after the FM detector and before the stereo decoder, marked with "80 mV". I think it would be possible to route the baseband signal from the test point to the input of a PC sound card for recording. But it's not required for a sound card to support the BW above 20 kHz, while all the decoding magic happens at the frequencies of up to 62 kHz. So I'm not sure the recorded file will be useful. Anyway, the stereo decoder in the receiver does not seem a very complex circuit. Perhaps you can get a clue by looking into the diagram.
I don't have an SDR dongle either. In my case, it just does not make sense. It's a big city and the reception is very poor. I can't hear anything on the AM bands no matter of the receiver. While the FM provides only a local content. So I'm not an SDR fan. I can only confirm the existence of the stereo signal in my area
-
It is a pity that you don’t have any SDR hardware available to make a direct recording of the signal. Even a very inexpensive dongle (10–20 USD) can do much more than just listening to radio – it allows spectrum analysis, interference studies from household equipment, debugging and aligning radio circuits across HF, VHF and UHF. I also live in the center of a large city with plenty of noise, but radio is still quite usable here: both HF and VHF broadcasting, amateur stations, satellites, and even frequency standard signals that can be used to calibrate frequency counters, generators, receivers or transceivers.
You are right that the stereo decoder circuitry itself is not extremely complex, but all the phase relationships and details matter a lot, so extracting the exact encoding parameters from schematics alone is unreliable – we really need real-world signals for verification. On your receiver’s diagram you indeed pointed out the correct place for recording the stereo composite: the KT7 80 mV test point. However, as you mentioned, capturing it with a standard sound card is problematic because of the input low-pass filter.
For composite recording an IQ signal is not required – it’s just a mono baseband with about 100 kHz bandwidth. In theory it could be digitized with a sound card at 384 kHz sampling rate (giving 192 kHz bandwidth), but in practice most sound cards have an input LPF around 20 kHz. A better option would be a capture card with 150–200 kHz effective bandwidth, or even a digital oscilloscope with sufficient memory depth. But I suspect you may not have such equipment either.
The most practical solution would be recording the IF with an SDR receiver. Perhaps you could borrow an RTL-SDR from a friend – these dongles are very popular among radio amateurs nowadays because they are inexpensive and versatile. They allow not only receiving nearly any band from longwave to UHF, but also testing and aligning radio equipment.
If you don’t have access to an SDR, you could still try recording the composite with a sound card, just in case. To compensate for the LPF effect, it would be very helpful if you could also record plain white noise from the same input – then I could attempt to restore the frequency response of your sound card and design a correction filter to recover as much of the high-frequency content as possible. Often the LPFs in sound cards are quite weak and do not completely suppress higher frequencies. Of course, this is a rather problematic and challenging approach to try to reconstruct the signal in such a way, but since there are no other available options to obtain an OIRT stereo composite sample, it might still be worth attempting.
It is especially regrettable that, from what you describe, your local station really does broadcast stereo in the OIRT system. According to the schematic, your receiver indeed detects OIRT stereo with the 31.25 kHz subcarrier, so there is little doubt about it. It is a pity that you don’t have the technical means to capture and record this signal. Broadcasts in this system are becoming quite exotic nowadays: much of the old transmitting equipment has aged and gone out of service, and modern replacements for OIRT seem to be very rare. In most cases stations still operating in this band transmit only in mono.
I am very interested in implementing an OIRT stereo encoder/decoder in DSP and experimenting with its pros and cons compared to the CCIR system. Unfortunately, in my region broadcasting in this band is only mono, so I have no chance to capture such signals myself.
Can you help with obtaining such OIRT stereo-composite recording for testing?
-
How many samples are needed? For me, the simplest option seems to capture the baseband signal with Rigol DS1054Z scope (check the specs. It's the only digitizing scope I have)
-
How many samples are needed? For me, the simplest option seems to capture the baseband signal with Rigol DS1054Z scope (check the specs. It's the only digitizing scope I have)
As many samples as possible :)
For example, with a Siglent scope that has 14 Mpts of memory, I can capture 14 seconds at a 1 MSa/s rate. To do this, it is necessary to manually disable Roll mode (which is enabled by default below 50 MSa/s) and to enable the 20 MHz low-pass filter to minimize aliasing. I think that 14 seconds should be sufficient to clearly hear the stereo.
There is, of course, a risk that aliases from frequencies above 500 kHz may affect the signal quality, but it is still worth trying.
-
Ok then. Take the receiver apart, blow the dust off and capture the signal. Will do it on the next week (hopefully). BTW for better anti-aliasing, I can arrange a simple passive LP filter (two res, two caps) with the cutoff frequency of about 100 kHz
-
Ok then. Take the receiver apart, blow the dust off and capture the signal. Will do it on the next week (hopefully). BTW for better anti-aliasing, I can arrange a simple passive LP filter (two res, two caps) with the cutoff frequency of about 100 kHz
No, it is better not to use an LPF at all than to add a simple 100 kHz one, since such a filter can significantly attenuate the amplitude of the encoded L–R signal, given that the OIRT stereo composite has a bandwidth of 100 kHz. If you do apply an LPF, it should have a much higher cutoff, around 300–400 kHz.
I am looking forward to the test sample — it will be interesting to see how the OIRT stereo composite actually looks.
-
I came across an article about polar modulation for stereo signals. It refers to a 31.25 kHz subcarrier, which suggests it is describing the OIRT stereo composite standard. However, what seems odd is that the article claims polar modulation assigns the left and right channels to the corresponding polarities of the subcarrier. This is unusual, since earlier publications indicated that the subcarrier carries the difference signal (L−R), in addition to the mono sum signal (L+R). Moreover, the detection circuit shown appears overly simplistic—essentially just two diodes.
What exactly is being described here? Is this an earlier revision of the standard, or simply a misunderstanding and the circuit itself is incorrect?
-
I haven't checked in detail yet but I think the described asymmetric AM does not fit within the +-15kHz bandwidth around the sub-carrier. And if you remove sub-carrier harmonics and their sidebands, I think you lose the intended asymmetry.
-
Yes, that’s exactly the confusing part. There are many references claiming that OIRT used some kind of “polar modulation”, but I suspect there is some kind of a misinterpretation. From what I can see, the system probably is essentially the same as CCIR stereo: the baseband carries L+R, and the 31.25 kHz subcarrier is DSB-modulated with L−R, only with the subcarrier left partially suppressed (at −14 dB) rather than fully suppressed as in CCIR.
The question then is: where do all the statements about direct diode decoding come from? The article above even includes a schematic with nothing more than two diode detectors and RC filters. This seems questionable — could such a simple circuit really decode a stereo signal in which the subcarrier is modulated with L−R?
That’s why I’m looking for an actual recording of an OIRT stereo broadcast. With a real composite signal it would be possible to analyze the spectrum and confirm what is really transmitted.
I did find a screenshot of an OIRT station transmitting in stereo — the composite signal spectrum is faintly visible there — but unfortunately that’s the only material I’ve managed to locate so far.
-
To radiolistener
-
Yes, that’s exactly the confusing part. There are many references claiming that OIRT used some kind of “polar modulation”, but I suspect there is some kind of a misinterpretation. From what I can see, the system probably is essentially the same as CCIR stereo: the baseband carries L+R, and the 31.25 kHz subcarrier is DSB-modulated with L−R, only with the subcarrier left partially suppressed (at −14 dB) rather than fully suppressed as in CCIR.
The question then is: where do all the statements about direct diode decoding come from? The article above even includes a schematic with nothing more than two diode detectors and RC filters. This seems questionable — could such a simple circuit really decode a stereo signal in which the subcarrier is modulated with L−R?
Doing some numerical experiments, it seems that that positive and negative rectification of a composite signal L + R + subcarrier DSB-modulated with L-R, can separate L and R. The separation was imperfect, though, and I did not get reasonable separation at all with DSB-RC. I also don't know if the experiment generalizes to arbitrary L and R signals, since the rectification is a non-linear operation, so additivity is not granted.
The asymmetric AM-modulation scheme described in the paper generates a composite signal containing a L+R component, the subcarrier DSB-modulated with L-R, and sidebands of subcarrier harmonics (but no subcarrier harmonics themselves). This signal can, of course, be separated to L and R by positive and negative rectification (and subsequent lowpass filtering), but due to the presence of sidebands of subcarrier harmormics, it has a much higher occupied bandwidth. After removing the sidebands of subcarrier harmormics with a lowpass, we would basically end up with a compositive signal L + R + subcarrier DSB-modulated with L-R (but not DSB-RC).
EDIT: Here's a simulation of the asymmetric AM with random L and R random signals, BW-limit to eliminate sidebands of sub-carrier harmonics, and demodulation by positive and negative rectification of the composite signal + lowpass (=> "ideal diode detector", but a much better lowpass). The demodulated signals are decomposed into L, R, DC components and residuals to assess the goodness of separation. With the random L and R the separation is not that bad. Std(resid) could be lower, I don't know yet where exactly it comes from.
pkg load signal
fs = 1e6
fsub = 31250
audiobw = 15000
SKIP = 20000
N = 100000 + SKIP
% generate bandwith-limited random L and R signals
L = 2 * rand(1,N) - 1;
R = 2 * rand(1,N) - 1;
[n,Wn,beta,ftype] = kaiserord([0.95*audiobw audiobw],[1 0],[0.0001 0.0001],fs);
hnoise = fir1(n,Wn,ftype,kaiser(n+1,beta),"noscale");
L = filter(hnoise,1,L);
R = filter(hnoise,1,R);
L /= max(abs(L));
R /= max(abs(R));
% L = cos(2*pi*2000*[0:N-1]/fs);
% R = cos(2*pi*3000*[0:N-1]/fs);
% subcarrier
subc = cos(2*pi*fsub*[0:N-1]/fs);
% asymmetric AM
composite = max(0,subc) .* (1 + L) + min(0,subc) .* (1 - R);
% limit BW of composite signal
[n,Wn,beta,ftype] = kaiserord([fsub+audiobw 2*fsub-audiobw],[1 0],[0.0001 0.0001],fs);
hcomp = fir1(n,Wn,ftype,kaiser(n+1,beta),"noscale");
composite = filter(hcomp,1,composite);
% lowpass after rectifier
[n,Wn,beta,ftype] = kaiserord([audiobw fsub-audiobw],[1 0],[0.0001 0.0001],fs);
hbb = fir1(n,Wn,ftype,kaiser(n+1,beta),"noscale");
LR1 = [L; R; ones(1,N)];
% compensate filter delays
delay = (length(hcomp)-1)/2 + (length(hbb)-1)/2
printf("Factor analysis of lowpass filtered composite signal:\n");
x = filter(hbb,1,composite);
lr1 = x(SKIP+1:end) * pinv(circshift(LR1,delay,2)(:,SKIP+1:end));
resid = x(SKIP+1:end) - lr1*circshift(LR1,delay,2)(:,SKIP+1:end);
printf(" mono: %9.6f * L + %9.6f * R + %9.6f + resid, std(resid) = %.6f\n", lr1, std(resid));
printf("Factor analysis of rectified + lowpass filtered composite signal:\n");
x = filter(hbb,1,max(0,composite));
lr1 = x(SKIP+1:end) * pinv(circshift(LR1,delay,2)(:,SKIP+1:end));
resid = x(SKIP+1:end) - lr1*circshift(LR1,delay,2)(:,SKIP+1:end);
printf(" positive: %9.6f * L + %9.6f * R + %9.6f + resid, std(resid) = %.6f\n", lr1, std(resid));
x = filter(hbb,1,min(0,composite));
lr1 = x(SKIP+1:end) * pinv(circshift(LR1,delay,2)(:,SKIP+1:end));
resid = x(SKIP+1:end) - lr1*circshift(LR1,delay,2)(:,SKIP+1:end);
printf(" negative: %9.6f * L + %9.6f * R + %9.6f + resid, std(resid) = %.6f\n", lr1, std(resid));
Example output:
fs = 1000000
fsub = 31250
audiobw = 15000
delay = 4012
Factor analysis of lowpass filtered composite signal:
mono: 0.317287 * L + 0.317287 * R + -0.000000 + resid, std(resid) = 0.000005
Factor analysis of rectified + lowpass filtered composite signal:
positive: 0.316714 * L + 0.000577 * R + 0.320410 + resid, std(resid) = 0.002654
negative: 0.000573 * L + 0.316710 * R + -0.320410 + resid, std(resid) = 0.002654
EDIT: My first conclusion:
- The described modulation scheme generates directly a composite signal (including the mono component), not just the modulated subcarrier.
- The unfiltered composite signal can be demodulated virtually exactly with two ideal diodes + LPF, but is has too much bandwidth
- After limiting the bandwith to 46.25 kHz, demodulation with two ideal diodes is still possible, but it is no longer perfect:
a) separation between L and R is no longer perfect, but still not too bad.
b) Besides L, R and DC components, the 0-15kHz spectrum contains some residual components(noise/distortion). Still not sure where it comes from; I guess IMD introduced by the rectification. The residual component is larger when L and R are two sine waves, and smaller when L and R are wideband signals.
- In the tube era, these flaws were likely considered acceptable.
- While it seems to be a working scheme (with some restrictions), I'm no sure if this is really OIRT. The article does not claim that either; it just talks about "Soviet Stereo". The modulated subcarrier amplitude is relatively high, and the subcarrier itself is not reduced and consumes almost 50% of the composite signal power. But exactly these amplitude ratios, as produced by the described modulator, seem to be necessary in order to make simple diode detection work.
-
After limiting the bandwith to 46.25 kHz, demodulation with two ideal diodes is still possible, but it is no longer perfect:
a) separation between L and R is no longer perfect, but still not too bad.
b) Besides L, R and DC components, the 0-15kHz spectrum contains some residual components(noise/distortion). Still not sure where it comes from; I guess IMD introduced by the rectification. The residual component is larger when L and R are two sine waves, and smaller when L and R are wideband signals.
Hm... so, if I understand correctly, it means that the stereo composite can indeed be demodulated using just two diodes, but the result will not be perfect in quality?
-
Hm... so, if I understand correctly, it means that the stereo composite can indeed be demodulated using just two diodes, but the result will not be perfect in quality?
To be more specific: It seems that this particular kind of stereo composite signal can be demodulated/decoded by ideal half-wave rectification + LPF (although the result is not perfect). How well the circuit can do the same job is a different question. This would require a circuit simulation. In particular, I cannot imagine an RC filter with a narrow 15...16.25 kHz transition band and high stop band attenuation at 16.25. But if you have speakers that cannot reproduce the residual ultrasound components in the output, and ears that cannot hear them, you may simply not care...
Btw, it's interesting that the described asymmetric AM method (plus limiting the bandwidth to 31.25+15 kHz) happens to generate a composite signal equivalent to
(L+R) / pi + subcarrier .* (1 + (L-R) / 2)
IOW, what we get is the mono signal (L+R)/2, scaled by factor 2/pi, plus the sub-carrier, DSB-modulated with (L-R)/2, 100% modulation. No reduced subcarrier. Although it seems to be a feasible method in principle, is this really OIRT Stereo? I don't think so.
EDIT:
... but the result will not be perfect in quality
Attached is an example of the demodulated R channel spectrum, when the original L and R are sine waves with 3kHz and 4kHz. SFDR is about 25dB, SINAD about 21dB. When L and R are wideband signals (random noise), I see a little bit better SINAD of roughly 29 dB.
EDIT: Added audio file of demodulated R channel, with 4 kHz tone.
-
Thanks to Njk, we obtained a test sample of an OIRT stereo composite signal, recorded with a Rigol oscilloscope connected at the test point between the FM demodulator output and the stereo decoder input of the Leningrad RP-015 receiver.
I applied a 50 kHz FIR low-pass filter for decimation and encoded the result into FLAC to reduce file size. The recording contains 24 seconds of music. The song is identified as "Prayer" by Tom Odell at 2:42: https://youtu.be/pTxeOOkTqYo?si=2aeP5-vgOQ81Fkii&t=162
The quality is fairly good, although some frequency distortions are audible. Applying a de-emphasis IIR filter significantly reduces them, but slight artifacts remain, most likely due to the limited linearity of the oscilloscope ADC. The 31.25 kHz pilot carrier is well visible.
The attached images show the spectrum of the original oscilloscope signal and the spectrum after applying the filter and decimation (file included).
Still not tried decoding it yet, only tried to apply a 50 µs IIR de-emphasis filter to the L+R component, which produced good results.
-
There can be some distortions because of multi-path reception. Very common phenomenon at my location. Tried to mitigate it by adjusting the antenna position but it's of dynamic nature and can happen at any time. BTW, another difference between OIRT and CCIR bands is that the former traditionally uses the horizontal polarization, so a receiving dipole is expected to be placed horizontally, like a TV antenna. While CCIR uses the vertical polarization. But for me, it makes no difference concerning the MPR. It's very annoying and happens with every receiver all the way
-
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.
-
The subcarrier level is about 100 mVpp. It can be measured when silence is broadcasting
-
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?
-
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)
[attach=1]
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:
[attach=2]
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:
[attach=3]
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:
[attach=4]
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
-
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):
[attach=1]
This can be expressed in Octave as:
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:
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:
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
-
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.
-
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?
-
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).
-
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*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.
-
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.
-
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 (https://www.eevblog.com/forum/rf-microwave/wfm-stereo-broadcasting-standards-oirt-vs-ccir/?action=dlattach;attach=2666767).
-
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.
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 (https://www.eevblog.com/forum/rf-microwave/wfm-stereo-broadcasting-standards-oirt-vs-ccir/?action=dlattach;attach=2666743;image) 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 (https://www.eevblog.com/forum/rf-microwave/wfm-stereo-broadcasting-standards-oirt-vs-ccir/?action=dlattach;attach=2666771) 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. ???
-
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 (https://www.eevblog.com/forum/rf-microwave/wfm-stereo-broadcasting-standards-oirt-vs-ccir/?action=dlattach;attach=2666743;image) 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 (https://www.eevblog.com/forum/rf-microwave/wfm-stereo-broadcasting-standards-oirt-vs-ccir/?action=dlattach;attach=2666771) 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).
-
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 ]
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/
-
Generated composite signal.
-
Stereo decoded from the generated composite, using the original subcarrier from the modulator.
-
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. ]
-
[ 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:
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 (https://liquidsdr.org/) library, but its less readable.
1)The MEX extension code
// 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:
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)
Cool, nice approach, I did it separately:
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?
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.
-
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
-
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.
-
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.
-
Can the gqrx stereo decoder demodulate Njk's recorderd composite w/o distortion?
-
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.
-
After a long pause, I returned to the OIRT multiplex demodulator. I finally made a working version for both CCIR and OIRT, and it seems to work ok :) for OIRT, the L/R channel separation is around 56 dB 8)
In GQRX it’s much worse — there are a few mistakes, and the PLL implementation in GNU Radio turns out to have insufficient step resolution, which leads to a fairly large phase error in carrier recovery. I had to rework the angular frequency integrator to use double precision, since float resolution wasn’t enough and it helps to keep pretty precise phase recovery.
The main issue with OIRT is the residual DC from the suppressed carrier after DSB-SC demodulation. In principle, it works fine on a synthetic signal even if I hardcode DC subtraction as 0.2, but in real stations the carrier level can drift slightly. So I added a DC estimator with an initial value of -0.2 and subtract the calculated DC level in real time. However, I don’t really like this approach — it may slightly cut low frequencies. :-\
I’m wondering if there’s a better way to reliably eliminate DC offset from carrier without affecting the low-frequency content?
-
Can the gqrx stereo decoder demodulate Njk's recorderd composite w/o distortion?
In my new demodulator, @Njk’s recorded OIRT samples do play, but with some distortion. At the same time, my own synthetic examples (generated by me) reproduce cleanly.
This strongly suggests the issue is with the nonlinearity of the oscilloscope ADC used for the recording. The voice sample sounds reasonably good overall, just with some added noise/interference. The music sample, however, shows clear distortion — the L−R component is noticeably affected.
I don’t have other real off-air OIRT recordings to compare yet, so for now I’m mostly testing with my own synthetic signals.
-
The main issue with OIRT is the residual DC from the suppressed carrier after DSB-SC demodulation. In principle, it works fine on a synthetic signal even if I hardcode DC subtraction as 0.2, but in real stations the carrier level can drift slightly. So I added a DC estimator with an initial value of -0.2 and subtract the calculated DC level in real time. However, I don’t really like this approach — it may slightly cut low frequencies. :-\
I’m wondering if there’s a better way to reliably eliminate DC offset from carrier without affecting the low-frequency content?
It's impossible to have a filter with a zero at DC without also affecting low frequencies near DC. However, when you talk about drift, it's not pure DC anyway, but rather low-frequency AC. Therefore, you actually want to remove these low frequencies as well. IMO, the only feasible solution is a high-pass filter. I think a cutoff of -0.1 dB at 20 Hz will be fine? The transmitted audio likely doesn't go down any further.
IIRs tend to become numerically unstable at such a large fs/fc ratio and they are not linear phase. FIR would need a large number of taps. So I'd tend to realize a complementary highpass as delay(signal) - lowpass(signal), where the lowpass is a simple 4th order moving average filter that can be realized with a CIC/RRS structure. With 1752-tap moving average (for each of the 4 stages), 3502 samples delay, and 48kSa/s, I get the attached magnitude response. That's a passband ripple of only ~0.02 dB and attenuation of 80 dB at 0.1 Hz (and linear phase, of course, so no phase distortion).
-
1752 taps feels quite excessive for a real-time SDR chain. I’d prefer to keep the demodulator lightweight — that kind of delay and buffer size is overkill for what is essentially slow drift / near-DC.
From my experiments, I ran into a different issue with the PLL. After switching from my own implementation to liquid-dsp nco_crcf, lock became much faster (~5x times), but I started getting audible spurs. After investigating, I found tones around the carrier at about -60 dBc at the PLL output.
I suspect this is due to the fixed loop parameters in liquid-dsp — there’s no control over damping, and the internal NCO/PLL seems to introduce phase modulation artifacts.
Also, OIRT demodulation turned out to be very sensitive to the phase detector. Using atan2(Q,I) - phase works fine, but other detectors break stereo demodulation, even though the PLL still locks. The reason isn’t fully clear yet for me.
So for now I’d prefer a simple/cheap DC removal and focus on getting a cleaner PLL/LO.
-
here is my attempt to implement PLL for Octave with using liquid-dsp nco_crcf in MEX file.
compile it with
mkoctfile --mex pll_process2_mex.cpp -lliquid
pll_process2_mex.cpp
// linux: mkoctfile --mex pll_process2_mex.cpp -lliquid
#include "mex.h"
#include <cmath>
#include <complex>
#include <liquid/liquid.h>
typedef std::complex<float> cf;
typedef struct {
nco_crcf nco;
float err;
float min_freq;
float 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.nco = nco_crcf_create(LIQUID_NCO);
nco_crcf_pll_set_bandwidth(pll.nco, loop_bw);
nco_crcf_set_frequency(pll.nco, 0.0f);
pll.err = 0.0;
pll.max_freq = max_freq;
pll.min_freq = min_freq;
}
static inline void pll_update(pll_t &pll, cf input) {
// phase detector
float phase = atan2f(input.imag(), input.real());
pll.err = mod_2pi( phase - nco_crcf_get_phase(pll.nco));
// update PLL
nco_crcf_pll_step(pll.nco, pll.err);
float freq = nco_crcf_get_frequency(pll.nco);
if (freq > pll.max_freq) freq = pll.max_freq;
if (freq < pll.min_freq) freq = pll.min_freq;
nco_crcf_set_frequency(pll.nco, freq);
nco_crcf_step(pll.nco);
}
//----------------------------------------------------------------------
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
if (nrhs != 4)
mexErrMsgTxt("Usage: [output, err] = pll_process2_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]);
// init
pll_t pll;
pll_init(pll, loop_bw, min_freq, max_freq);
for (mwSize n = 0; n < N; n++) {
cf lo;
nco_crcf_cexpf(pll.nco, &lo);
y_re[n] = lo.real();
y_im[n] = lo.imag();
err_re[n] = pll.err;
pll_update(pll, cf(x_re[n], x_im[n]));
}
}
it requires to install liquid-dsp library: https://github.com/jgaeddert/liquid-dsp.git
Here is Octave test script:
pkg load signal;
fs = 384000; % sample rate
f0 = 19000; % CCIR pilot
N = fs/8; % 0.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));
% --- FIR LPF and complex BPF kernel calculator ---
function ntaps = compute_ntaps(fs, transition_width, max_atten)
ntaps = (max_atten * fs) / (22.0 * transition_width);
ntaps = floor(ntaps);
% --- ensure odd length (linear-phase FIR requirement)
if mod(ntaps, 2) == 0
ntaps = ntaps + 1;
end
end
function taps = low_pass(gain, fs, fc, tw)
% construct the truncated ideal impulse response
% [sin(x)/x for the low pass case]
ntaps = compute_ntaps(fs, tw, 53);
w = hamming(ntaps).';
M = floor((ntaps - 1)/2);
n = -M:M;
wc = 2*pi*fc/fs;
taps = zeros(1, ntaps);
idx = (n == 0);
taps(idx) = (wc/pi) * w(idx);
idx = ~idx;
taps(idx) = (sin(n(idx)*wc) ./ (n(idx)*pi)) .* w(idx);
% normalize
fmax = taps(M+1) + 2*sum(taps(M+2:end));
gain = gain / fmax;
taps = taps * gain;
end
function taps = complex_band_pass(gain, fs, f1, f2, tw)
% low-pass prototype
lptaps = low_pass(gain, fs, (f2-f1)/2, tw);
ntaps = length(lptaps);
n = 0:ntaps-1;
fc = (f1 + f2)/2;
w0 = 2*pi*fc/fs;
% modulation
taps = lptaps .* exp(1j * w0 * (n - (ntaps-1)/2));
end
% --- initialize PLL ---
pll_bw = 10;
pll_minf = f0-50;
pll_maxf = f0+50;
pilotc = complex(pilot*2.0, 0);
bpf = complex_band_pass(1.0, fs, pll_minf, pll_maxf, 5000.0);
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);
%plotFreqz(bpf, 1, fs, -fs/2, fs/2, -150, 0, 2^20, true);
pilotc = filter(bpf, 1, pilotc);
fprintf('bpf done\n');
[pll_lo, err] = pll_process2_mex(pilotc, pll_bw * 2.0*pi/fs, pll_minf * 2.0*pi/fs, pll_maxf * 2.0*pi/fs);
fprintf('pll done\n');
% --- align with BPF delay ---
pilotc(end+1:end+bpf_delay) = zeros(1, bpf_delay);
pll_lo(end+1:end+bpf_delay) = zeros(1, bpf_delay);
err(end+1:end+bpf_delay) = zeros(1, bpf_delay);
pilot = [zeros(bpf_delay, 1); pilot];
% --- plot FFT Spectrum with Kaiser window ---
function plot_fft_kaiser(x, fs)
Nfft = length(x);
beta = 8.6; % ~80 dB sidelobes
w = kaiser(Nfft, beta);
xw = x(:) .* w; % apply window
X = fftshift(fft(xw));
% normalize (avoid log(0))
mag = abs(X);
mag = mag / max(mag + 1e-12);
mag_db = 20 * log10(mag + 1e-12);
f = (-Nfft/2:Nfft/2-1) * (fs / Nfft);
figure;
plot(f, mag_db);
grid on; grid minor;
set(gca, 'gridalpha', 0.2); set(gca, 'minorgridalpha', 0.05);
xlabel('Frequency (Hz)');
ylabel('Magnitude (dB)');
xlim([-fs/2 fs/2]);
ylim([-160 5]);
end
% skip unstable part
pilot_fft = pll_lo(1000:end);
pot = 2 ^ floor(log2(length(pilot_fft)));
pilot_fft = pilot_fft(1:pot);
plot_fft_kaiser(pilot_fft, fs);
% --- plot Waveform ---
N = length(pilotc);
t = 1:N;
figure;
hold on;
plot(t, real(pilotc), 'DisplayName', 'pilotc');
plot(t, real(pll_lo), 'DisplayName', 'pll_lo');
plot(t, err, 'DisplayName', 'err');
grid on; grid minor;
set (gca, 'gridalpha', 0.2); set(gca, 'minorgridalpha', 0.05);
xlim([min(t) max(t)]);
legend();
Is there a way to eliminate these spurs?
-
Is there a way to eliminate these spurs?
Good question where they come from. Can you plot pll.err over time? Zoom in and check if it ripples quicky, like one sample up, one sample down, and so on? Does it help if you dcrease loop BW by factor 10 or maybe even 100?
-
NCO spurs due to phase truncation, maybe. What kind of precision do you get from the LiquidDSP oscillator implementation?
-
Is there a way to eliminate these spurs?
Good question where they come from. Can you plot pll.err over time? Zoom in and check if it ripples quicky, like one sample up, one sample down, and so on? Does it help if you dcrease loop BW by factor 10 or maybe even 100?
I checked that. Reducing loop BW actually makes the ripple around the spurs more pronounced, but the spurs themselves don’t go away.
There is some ripple visible on pll.err, but it’s roughly at the same level as with my manual PLL implementation — and that version doesn’t produce any spurs.
So it doesn’t look like the issue is simply error signal noise or loop bandwidth. It seems more like something intrinsic to the NCO/PLL implementation rather than the loop bandwidth.
-
NCO spurs due to phase truncation, maybe. What kind of precision do you get from the LiquidDSP oscillator implementation?
From what I see in the liquid-dsp nco_crcf implementation, it uses a 32-bit integer phase accumulator, i.e. effectively a 2^32 phase resolution.
// type == LIQUID_NCO, LIQUID_VCO_INTERP
uint32_t theta; // 32-bit phase [radians]
uint32_t d_theta; // 32-bit frequency [radians/sample]
The output is then generated via a lookup table (or interpolated LUT depending on configuration).
Here is the source for pll_set and pll_step from liquid-dsp source: liquid-dsp/src/nco/src/nco.proto.c:
// set pll bandwidth
int NCO(_pll_set_bandwidth)(NCO() _q,
T _bw)
{
// ...
_q->alpha = _bw; // frequency proportion
_q->beta = SQRT(_q->alpha); // phase proportion
// ...
}
// advance pll phase
// _q : nco object
// _dphi : phase error
int NCO(_pll_step)(NCO() _q,
T _dphi)
{
// ...
// increase frequency proportional to error
NCO(_adjust_frequency)(_q, _dphi*_q->alpha);
// TODO: ensure frequency doesn't run away
// increase phase proportional to error
NCO(_adjust_phase)(_q, _dphi*_q->beta);
// constrain frequency
//NCO(_constrain_frequency)(_q);
// ...
}
// adjust frequency of nco object
int NCO(_adjust_frequency)(NCO() _q,
T _df)
{
// ...
_q->d_theta += NCO(_constrain)(_df);
// ...
}
// adjust phase of nco object, constraining phase
int NCO(_adjust_phase)(NCO() _q,
T _dphi)
{
// ...
_q->theta += NCO(_constrain)(_dphi);
// ...
}
-
KE5FX is right. I get the same phase truncation spurs when I calculate DDS (free-running, no PLL) manually, from only 1024 wavetable entries. 1024 entries with phase truncation are simply not enough if you want a higher SFDR. You need either more entries (larger NCO_STATIC_LUT_NBITS), or you interpolate, or both. Type LIQUID_VCO_INTERP seems to interpolate.
fs = 384000
fc = 19000
wt = single(exp(2i*pi*[0:1023]/1024));
fcw = round(fc/fs*2^32)
phi = floor([0:fs-1]*fcw/2^22);
sig = wt(1+floor(mod(phi,1024)));
w = flattopwin(length(sig),"periodic")';
plot(([0:fs-1]-fs/2)*0.001,20*log10(abs(fftshift(fft(sig.*w)/sum(w))))); grid on
Btw, I deleted my previous post. There was an off-by-one error in the calculation (and consequently, also in the conclusion).
EDIT: With a 64k wavetable and truncation, I get them down to -96 dBc, see figure 2.
EDIT: Seems that 1024 entries with linear interpolation are even more effective than 64k with truncation, see figure 3.
fs = 384000
fc = 19000
wt = single(exp(2i*pi*[0:1024]/1024));
phi = mod([0:fs-1]*fc/fs*1024,1024);
sig = interp1([0:1024],wt,phi);
w = flattopwin(length(sig),"periodic")';
plot(([0:fs-1]-fs/2)*0.001,20*log10(abs(fftshift(fft(sig.*w)/sum(w))))); grid on
-
That’s the interesting part — in my manual PLL implementation this code:
lo = cexpf(I * pll.phase);
does not produce visible spurs.
However, when I use the same oscillator approach with the liquid-dsp PLL phase state instead of nco_crcf_cexpf(pll.nco, &lo), spurs still appear, but their position shifts closer to the carrier and with less magnitude.
-
That’s the interesting part — in my manual PLL implementation this code:
lo = cexpf(I * pll.phase);does not produce visible spurs.
// get phase [radians]
T NCO(_get_phase)(NCO() _q)
{
if (_q->type == LIQUID_VCO_DIRECT) {
return liquid_error(LIQUID_EICONFIG,"error: nco_get_phase(), "
"cannot be used with object type == LIQUID_VCO_DIRECT");
}
return TIL(2)*TFL(M_PI)*(T)_q->theta / (T)(1LLU<<32);
}
Note that _q->theta is the phase at full 32-bit precision of the phase accumulaor, not yet truncated to its upper 10 bits, which act as index into the wavetable.
[ phase range [0...2*pi) maps to theta range [0...2^32) ]
Your cexpf(I * pll.phase) keeps ~24-bit precision (single precision float) for the phase, unlike the 10-bit table lookup.
Did you already try LIQUID_VCO_INTERP?
-
Here is the current PLL behavior (my implementation) on a real OIRT signal previously posted by Njk in this thread.
Does this look like normal PLL operation, or is something wrong?
The large phase error fluctuations are caused by the presence of DSB from L-R components after applying the BPF to composite.
Parameters:
fs = 96 kHz
PLL loop bandwidth = 5 Hz
BPF 183 taps:
f0 = 31250 Hz
cutoff = ±50 Hz
transition width = 5000 Hz
stop-band attenuation = 60 dB
If the BPF transition width is reduced to 100 Hz, the phase error fluctuations disappear and the PLL produces a nice smooth curve. However, the BPF length then grows to about 9065 taps, which is impractical. Increasing the BPF from 183 to 9065 taps does not noticeably improve stereo separation — the only visible difference is that the lock metric becomes cleanly equal to 1.
The current demux implementation provides about 60 dB L/R channel separation for OIRT and about 100 dB for CCIR. However, I suspect the limitation may come from the de-emphasis response curve, due to slight phase/amplitude mismatch between the pre-emphasis and de-emphasis IIR filters. This still needs to be tested.
-
Do you still use atan phase detector? Then I guess this might be the problem. It needs a clean, pre-filtered pilot. A mulitiplying detector (mixer) should not need it -- the loop filter does the filtering.
It can't lock half vco frequency to pilot, though, so you may need an NCO that generates pilot sin (for pd) and cos (for stereo detection) and subcarrier cos waves simultaneously.
EDIT: The last sentence is, of course, not relevant for OIRT, where we don't have a separate pilot tone with half subcarrier frequency.
-
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));
-
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.
Thinking more deeply about this, it seems CCIR has a solid ±4 kHz gap around the pilot tone, which is PLL-friendly. OIRT is much tighter: assuming audio goes down to 30 Hz, you only have a ±30 Hz gap on each side of the subcarrier before the L-R audio subbands start. Consequently, this requires narrow pre-filtering in front of the PLL and/or a tight PLL loop bandwidth. This does not pose a problem for slow frequency drift, but it reduces the system's tolerance for following FM-induced subcarrier jitter, and locking at larger initial frequency offsets.
-
Yes, I see the same thing — in CCIR the carrier is recovered much more cleanly because of the wider spectral gap around the pilot.
But the interesting part is that with the atan2-based detector, OIRT still works even with a very wide BPF transition width of about 5 kHz, while other phase detectors fail in that case. So apparently the atan2 detector is much less sensitive to the DSB components present around the carrier.
The main issue is that a moderate reduction of the BPF transition width does not produce any noticeable improvement. Significant improvement only starts to appear when the transition width is narrowed to around 100 Hz or less. However, this increases the filter length to more than 9000 taps, which is extremely expensive from a performance standpoint.
From what I can see, this would suggest that a very narrow filter is required. But at the same time the system appears to work reasonably well even with a 5 kHz transition width when using the atan2-based phase detector. So at first glance this does not seem like a problem that justifies such a massive increase in computational cost.
Are there any alternative approaches to improve this behavior without requiring such an extremely long BPF?
-
Which composite signal? Can you repost the audio file, or refer to the msg where it was already atached?
-
Which composite signal? Can you repost the audio file, or refer to the msg where it was already atached?
I can't attach the original samples in good quality due to the attachment size limit.
Even with FLAC 24-bit compression, the file is 14 MB, while the limit here is just 8 MB.
Here is a link to the file-sharing host. They say the file should be available for 10 days:
- test-oirt-1-lpf96k.wav: https://gofile.io/d/DKwfAp
- test-oirt-2-lpf96k.wav: https://gofile.io/d/1VsR6s
Unfortunately, the quality is not very good since the signal was captured using an oscilloscope, but these are currently the only real off-air OIRT composite samples I have available.
The second sample (music) demuxes to stereo, but with noticeable frequency distortions. The first sample (voice) demuxes quite well overall, although there is some slight background noise.
-
Here is an additional my synthetic sample that I generated in Octave.
But note, during generation, a pre-emphasis IIR filter H(z) = 1 + k - k*z^-1; where k = Fs * tau; (as used in CCIR systems) was applied instead of the frequency-shaping formula specified in the OIRT standard. Therefore, this sample should not be used as a reference, but it works and demultiplexes in stereo.
-
Yes, I see the same thing — in CCIR the carrier is recovered much more cleanly because of the wider spectral gap around the pilot.
But the interesting part is that with the atan2-based detector, OIRT still works even with a very wide BPF transition width of about 5 kHz, while other phase detectors fail in that case. So apparently the atan2 detector is much less sensitive to the DSB components present around the carrier.
The main issue is that a moderate reduction of the BPF transition width does not produce any noticeable improvement. Significant improvement only starts to appear when the transition width is narrowed to around 100 Hz or less. However, this increases the filter length to more than 9000 taps, which is extremely expensive from a performance standpoint.
From what I can see, this would suggest that a very narrow filter is required. But at the same time the system appears to work reasonably well even with a 5 kHz transition width when using the atan2-based phase detector. So at first glance this does not seem like a problem that justifies such a massive increase in computational cost.
Just a few thoughts on this. Trying to approach it analytically, it seems that an atan2 phase detector is well-suited for an AM-modulated analytic signal
(1 + modulating_sig(t)) · e^(i·(ω·t + φ)) [ where ω is the carrier frequency ]
because atan2 is amplitude-independent. As long as modulating_sig(t) > -1 is granted, the following applies:
angle( (1 + modulating_sig(t)) · e^(i·(ω·t + φ)) · e^(-i·ω·t) ) = φ [ independent of modulating_sig :) ]
However, if modulating_sig(t) > -1 is not met, the detector suffers from 180° phase jumps/inversions (i.e. discontinuities) in the phase_error(t) signal.
Thus, for the atan2 phase detector to work properly with an AM signal, three conditions are IMO required:
1. The signal must be analytic (no negative frequencies).
2. Extraneous components like L+R must be filtered out.
3. Over-modulation must be avoided.
Conditions (1) and (2) are achieved if your complex bandpass filter has sufficient attenuation for frequencies < 15 kHz. However, condition (3) does not automatically follow. Once the subcarrier is reduced, the phase detector suffers from over-modulation (wrt. the reduced subcarrier) unless L−R happens to be weak enough.
And now the ominous transfer function K(f) from the standard comes into play. It specifies how the carrier and the sidebands near the carrier are reduced. Consequently, applying the inverse transfer function should undo this reduction and fully restore the original AM signal — complete with its full carrier and free of over-modulation.
While K(f) is defined as a shelved high-pass filter in the baseband, which is applied before modulation, it could equivalently be frequency-translated to the subcarrier frequency and applied after regular AM modulation. The same applies to the inverse of K(f) — it can also be frequency-shifted to the subcarrier frequency and applied before demodulation.
Basically, I think that any bandpass filter whose peak-aligned frequency response stays below the curve of the frequency-translated inverse of K(f) should achieve the goal of eliminating over-modulation. This leads us to a maximum bandwidth of ~324 Hz for the bandpass filter, assuming a complex 1st-order shape. Then, the attenuation at 15kHz becomes ~40dB, which I think is sufficient as well.
Are there any alternative approaches to improve this behavior without requiring such an extremely long BPF?
What about a 1st order complex IIR bandpass?
At resonance, the phase is exactly 0°. Sure, a frequency offset of the subcarrier introduces a static phase shift, and any phase offset of the recovered subcarrier translates into a cosine amplitude error in the demodulated L-R path, which in turn limits channel separation. However, the impact is likely acceptable. Even with a 5° phase offset, the channel separation remains at over 54 dB. With a bandpass bandwidth of ~300 Hz, a 5-degree phase budget maps to a frequency tolerance of roughly ±450 to ±500 ppm.
-
Are there any alternative approaches to improve this behavior without requiring such an extremely long BPF?
Yet another idea:
At (say) 10 Hz loop bandwidth, updating the PLL 384,000 times per second is overkill.
Doing 100–200 updates per second should be enough to sample the error signal properly.
Here is the proposed approach:
1. Multiply the (real) input stream with the conjugate of the complex NCO (down-mix).
2. Divide the down-mixed stream into blocks of 2048 samples each.
3. Apply a 2048-tap (periodic) Hann window to each block and sum the complex samples up.
4. Calculate the phase error = angle(sum) at the end of each block to update the PLL.
Assuming a loop bandwidth of 10 Hz, the additional intra-loop delay of approximately 2.7 ms (1024 samples) degrades the phase margin by ~9°. This reduction is IMO acceptable and should not induce loop oscillations. Although alternative window functions or FIR kernels could be used in step 3, I think that a Hann window is a good trade-off. It has a 60 dB/decade asymptotic roll-off, ensuring more than 100 dB attenuation at <= 15 kHz, despite a narrower equivalent noise bandwidth compared to other windows (except rectangular which is not an option anyway). Moreover, with 2048 taps, the resulting -3 dB bandwidth is ~270 Hz. Since this is narrower than the inverse of K(f), I think it will suffice to "over-undo" K(f) in order to eliminate over-modulation.
EDIT: Unlike an out-of-loop pre-filter, this filter is embedded inside the loop.
Consequently, the VCO remains in-phase with the input signal and suffers no filter time delay.
EDIT: Naturally, the PLL loop filter parameters must be adjusted to account for the lower fs/2048 update rate (~187.5 Hz) to maintain the target 10 Hz loop bandwidth.
-
Here is an additional my synthetic sample that I generated in Octave.
My own PLL attempt can lock to your sample virtually perfectly. Residual ripple in the phase detector output is less than a millidegree (before loop filter).
The other two recorded samples also do lock-in and stay locked. I don't see phase jumps and inversions.
But ripple in phase detector output is up to 10-20° peak :(
I've no idea yet what specific "anomaly" of the signal is the root cause :-//
I have not yet tried to demodulate.
-
I found a synthetic recording which appears to have been generated by Stereo Tool, and it seems to produce an incorrect OIRT multiplex. I do not see the expected polar channel separation in the waveform.
However, the file is interesting for another reason. My decoder correctly demodulates the initial stereo fragments that are present only in the left or only in the right channel, i.e. when the L-R signal appears. But at the same time, the PLL lock detector suddenly collapses because the phase detector output starts producing extremely large swings, almost reaching min/max excursions.
I do not fully understand what could produce such huge phase error spikes. The file is clearly generated incorrectly in some way, but it is not obvious what exactly is wrong with it that causes the phase detector to behave so violently.
What is even more interesting is that if the BPF transition width is narrowed to about 100 Hz, the PLL locks cleanly and the phase error spikes disappear completely. So the issue is clearly related to the nearby SSB sidebands somehow severely disturbing the phase estimate.
It is interesting to understand why these sidebands are able to corrupt the detected phase so strongly in this particular signal.
As for this strange file, a look at the waveform clearly shows that there are issues with the carrier amplitude, and the characteristic polar channel separation is missing. That polar separation is one of the key features of the OIRT multiplex, so its absence strongly suggests that the signal was generated incorrectly. But the main question here is what exactly is happening inside the phase detector itself (mix-down + atan2). That is the part which currently looks the most puzzling.
You mentioned earlier that over-modulation could potentially cause this kind of behavior. Could this perhaps be exactly such a case?
As far as I understood from the Stereo Tool author’s comments, he attempted to implement the frequency-shaping formula described in the standard. Perhaps this is exactly where the issue originates. What do you think?
-
I found a synthetic recording which appears to have been generated by Stereo Tool, and it seems to produce an incorrect OIRT multiplex. I do not see the expected polar channel separation in the waveform.
However, the file is interesting for another reason. My decoder correctly demodulates the initial stereo fragments that are present only in the left or only in the right channel, i.e. when the L-R signal appears. But at the same time, the PLL lock detector suddenly collapses because the phase detector output starts producing extremely large swings, almost reaching min/max excursions.
I do not fully understand what could produce such huge phase error spikes. The file is clearly generated incorrectly in some way, but it is not obvious what exactly is wrong with it that causes the phase detector to behave so violently.
What is even more interesting is that if the BPF transition width is narrowed to about 100 Hz, the PLL locks cleanly and the phase error spikes disappear completely. So the issue is clearly related to the nearby SSB sidebands somehow severely disturbing the phase estimate.
It is interesting to understand why these sidebands are able to corrupt the detected phase so strongly in this particular signal.
These are clearly phase inversions that are still left over after filtering.
Are you sure that the K(f) filter, as specified in the standard, was correctly applied to the L-R signal prior to subcarrier modulation? Because K(f) is a shelved highpass filter that attenuates low frequencies by up to 14 dB, it naturally attenuates the close-in sidebands near the subcarrier correspondingly.
Once your bandpass filter becomes narrower than the lowest audio frequency in the recording, you end up eliminating the sidebands completely, which explains why this approach works. However, with proper transmitter-side K(f) filtering, I think that a ~300 Hz bandpass should generally suffice, independent of the lowest audio frequency. BTW, the bandpass response should be bell-shaped rather than flat-topped, as the latter is counterproductive for interaction with the K(f) frequency response.
EDIT: Deletet potentially wrong statement - need to re-think.
EDIT: Regarding my point that "a 300 Hz bandpass should generally suffice": it is not merely a matter of bandwidth, but also of the filter's shape. The condition for atan2 phase detector stability is | 5 * K(f) * HLP(f) | <= 1, for any frequency f, where HLP represents the frequency response of the lowpass prototype of the filter in front of the phase detector. When using a Hann window (what I do), the minimum number of taps required is 0.003974 * Fs (-> 763 @192 kSa/s). Prerequisite is, of course, that the transmitter does K(f) filtering.
The worst-case test signal is certainly one where R = -L applies, at 100% level, but it is not possible to specify a particular worst-case audio frequency, since for an arbitrary HLP(f), the condition | 5 * K(f) * HLP(f) | could become >= 1 at any frequency.
Here's what I assume on the TX side:
L = ...
R = ...
% TODO: pre-emphasis - not relevant for this consideration
% K(f) shelved highpass, reduces the subcarrier and close-in sidebands
tau = 4.6 / (2000 * pi);
[b,a] = bilinear([tau 1], [tau 5], 1/Fs);
composite = 0.8*(L+R)/2 + filter(b,a,1+0.8*(L-R)/2) .* cos(2*pi*fsub*t);
-
After re-reading the standard more carefully (see the attached screenshot), I noticed a few details that may be relevant here.
According to sections 2.1.2.1 and 2.1.2.2, pre-emphasis is explicitly applied to both M=(L+R)/2 and S=(L-R)/2 before the subcarrier modulation stage. So I am not fully convinced that pre-emphasis can be ignored for this discussion, because it changes the spectrum of S before the shaping process takes place.
The other thing that still looks somewhat ambiguous to me is the interpretation of K(f). The standard states that the spectrum of the amplitude-modulated subcarrier shall be transformed according to K(f), but it does not explicitly describe the implementation method itself. It specifies the desired spectral result, but not exactly how to achieve it.
That leaves some room for interpretation:
- apply K(f) directly to S before modulation,
- apply shaping after AM modulation,
- or use some equivalent implementation.
These approaches may produce similar spectra while not necessarily producing identical phase behavior, which could potentially matter for the phase detector.
Regarding the strange file itself: it was generated by Stereo Tool, but the OIRT implementation there appears to have been only an early experimental attempt. The author spent quite some time trying to understand what K(f) actually means and eventually implemented it as an IIR approximation of some sort. However, I am not sure it was ever verified against the standard, and I am not even sure whether the version used to generate this file contained that implementation at all.
As far as I can tell, no further work or discussion happened afterwards, so it is entirely possible that the generated multiplex is not fully OIRT compliant.
Discussion thread about OIRT in Stereo Tool: https://forums.stereotool.com/viewtopic.php?t=6233
There may be a typo: 4.6 was likely intended to be 6.4. Would this be the correct form?
L = ...
R = ...
M = (L+R)/2; % L+R
S = (L-R)/2; % L-R
% ... pre-emphasis for M and for S
% K(f) shelved highpass, reduces the subcarrier and close-in sidebands
tau = 6.4 / (2*pi*1000);
[b,a] = bilinear([tau 1], [tau 5], 1/Fs);
subcf = 31250; % subcarrier frequency 31.25 kHz
composite = 0.8*M + filter(b,a,1+0.8*S) .* cos(2*pi*subcf*t);
-
Regarding K(f), I have a question about its practical impact on the system behavior.
It seems that applying K(f) increases the relative level of the close-in SSB components with respect to the carrier, and this may potentially interfere with the atan2-based phase detector.
I am considering introducing an inverse K(f) filter for spectral correction, since it also affects the reconstruction of the L-R channel.
However, this raises a couple of structural questions about the processing chain:
My current pipeline is:
IF → BPF → mixdown → atan2 → NCO output
and in parallel:
IF → mixdown → DSB demod → LPF → L-R extraction
Where would the inverse K(f) be correctly applied in this structure?
Should it be applied before both branches, i.e. something like:
IF → mixdown → K⁻¹(f) → mixup → (then both PLL and L-R processing paths)?
Or is there a more correct placement?
Since K(f) has a non-trivial frequency-dependent phase response, it affects not only magnitude but also the phase structure of the L-R component.
Does this imply that the inverse filter must reproduce both magnitude and phase characteristics of K(f) very accurately in order to preserve correct stereo separation? Otherwise, could this introduce additional separation errors between channels?
Am I understanding this correctly?
-
After re-reading the standard more carefully (see the attached screenshot), I noticed a few details that may be relevant here.
According to sections 2.1.2.1 and 2.1.2.2, pre-emphasis is explicitly applied to both M=(L+R)/2 and S=(L-R)/2 before the subcarrier modulation stage. So I am not fully convinced that pre-emphasis can be ignored for this discussion, because it changes the spectrum of S before the shaping process takes place.
I see pre-emphasis as relevant only with regard to the maximum levels of M, S and composite signal. If S becomes too strong after pre-emphasis, you would over-modulate the subcarrier. Either you ensure proper levels with a proper static gain for L and R, or with a dynamic compressor. Radio stations certainly do the latter. Otherwise, I don't see how pre-emphasis matters for the consideration.
The other thing that still looks somewhat ambiguous to me is the interpretation of K(f). The standard states that the spectrum of the amplitude-modulated subcarrier shall be transformed according to K(f), but it does not explicitly describe the implementation method itself. It specifies the desired spectral result, but not exactly how to achieve it.
That leaves some room for interpretation:
- apply K(f) directly to S before modulation,
- apply shaping after AM modulation,
- or use some equivalent implementation.
The specified transfer function is a baseband TF (before modulation). It can be realized with a simple passive voltage divider where the upper arm is R1 || C, and the lower arm is R2, where R1 = 4 * R2, and tau = R1 * C.
Any equivalent implementation is certainly allowed. But an exact post-modulation implementation using an analog shelved notch filter is actually impossible, because the symmetry of analog (or digital IIR) bandpass and notch filters is geometric rather than arithmetic. Good question how much asymmetry and thereforel AM -> PM conversion an analog 2nd order shelved notch filter would introduce.
These approaches may produce similar spectra while not necessarily producing identical phase behavior, which could potentially matter for the phase detector.
Crucial for the instantaneous phase is the conjugate complex symmetry between the upper and lower modulation sidebands, as any asymmetry introduces phase modulation in addition to AM.
There may be a typo: 4.6 was likely intended to be 6.4.
Thanks for pointing out the typo.
-
Regarding K(f), I have a question about its practical impact on the system behavior.
It seems that applying K(f) increases the relative level of the close-in SSB components with respect to the carrier, and this may potentially interfere with the atan2-based phase detector.
n of the L-R channel.
However, this characteristic is precisely what enables the atan2 phase detector to operate reliably at all, without requiring an extremely narrow pre-filter.
However, this raises a couple of structural questions about the processing chain:
My current pipeline is:
IF → BPF → mixdown → atan2 → NCO output
and in parallel:
IF → mixdown → DSB demod → LPF → L-R extraction
...
Where would the inverse K(f) be correctly applied in this structure?
I would suggest the following structure:
composite(real) → M'
composite(real) → mixer(real x complex) → lowpass → [ optionally down-sample ] → atan2 → PLL loop filter → VCO(complex) → back to mixer
composite(real) → mixer(real x complex) → extract real part → S'
M' → 15 kHz lowpass → de-emphasis → M
S' → K-1(f) → 15 kHz lowpass → de-emphasis → S
The ordering of cascaded filters is irrelevant in linear system theory; however, structural rearrangement is often driven by practical implementation constraints. So feel free to re-arrange. Furthermore, because K(s) · K⁻¹(s) = 1, introducing a compensation delay in the M path is unnecessary, provided that K⁻¹(s) is realized as an analog or digital IIR filter.
In the flowchart above, I replaced the band-pass filter in front of the loop with a low-pass filter inside the PLL loop. For a loop bandwidth of less than 10 Hz, the loop should still remain stable despite the phase margin reduction introduced by the low-pass filter.
The optional down-sampling step reduces the computational overhead of the atan2 operation and PLL updates by eliminating the need for per-sample execution. Additionally, if the low-pass stage is a FIR filter, the computational complexity is significantly minimized when the filter calculation is combined with down-sampling. This can be implemented as a polyphase decimator, or, if the down-sampling factor equals the filter length, it even becomes a highly efficient block filter.
Particularly in the block filter configuration, the FIR approach is significantly more efficient than an IIR filter. Nevertheless, I am still investigating whether a 3rd-order critically damped IIR filter could be used for the non-down-sampled case.
Since K(f) has a non-trivial frequency-dependent phase response, it affects not only magnitude but also the phase structure of the L-R component.
Does this imply that the inverse filter must reproduce both magnitude and phase characteristics of K(f) very accurately in order to preserve correct stereo separation? Otherwise, could this introduce additional separation errors between channels?
Yes, I would say so. In the analog domain, you can undo it exactly using a simple R1, R2, C filter and a 5x gain stage. In the digital domain, you need to check how well the magnitude and phase responses of the bilinear transform of K⁻¹(s) match those of the continuous-time K⁻¹(s). Maybe K⁻¹(z) needs to be tweaked a bit. [ Well, I said 'exact,' but in reality, analog also suffer from component tolerances. ]
EDIT:
Nevertheless, I am still investigating whether a 3rd-order critically damped IIR filter could be used for the non-down-sampled case.
The following IIR filter also seems to work ffine for me.
Tau is the same as for K(f). We need 3rd order, in order to get > 100 dB for 15kHz mono components.
% phase detector pre-filter
[b3,a3] = bilinear(1, [tau^3, 3*(tau^2), 3*tau, 1], 1/Fs);
-
I would suggest the following structure:
composite(real) → M'
composite(real) → mixer(real x complex) → lowpass → [ optionally down-sample ] → atan2 → PLL loop filter → VCO(complex) → back to mixer
composite(real) → mixer(real x complex) → extract real part → S'
M' → 15 kHz lowpass → de-emphasis → M
S' → K-1(f) → 15 kHz lowpass → de-emphasis → S
The ordering of cascaded filters is irrelevant in linear system theory; however, structural rearrangement is often driven by practical implementation constraints. So feel free to re-arrange. Furthermore, because K(s) · K⁻¹(s) = 1, introducing a compensation delay in the M path is unnecessary, provided that K⁻¹(s) is realized as an analog or digital IIR filter.
In the flowchart above, I replaced the band-pass filter in front of the loop with a low-pass filter inside the PLL loop. For a loop bandwidth of less than 10 Hz, the loop should still remain stable despite the phase margin reduction introduced by the low-pass filter.
One point still confuses me a bit.
In your proposed structure, K⁻¹(f) is applied only in the S' path for L-R reconstruction, while the PLL path still operates on the original K(f)-shaped signal:
composite → mixdown → LPF → atan2 → PLL
Wouldn't this leave the PLL phase detector exposed to the same altered carrier/sideband ratio introduced by K(f)?
My concern is that without K⁻¹(f) before the PLL path, the close-in sidebands around the subcarrier would still remain relatively elevated with respect to the carrier. In that case it seems that the PLL branch would still require a very narrow LPF, about ~100-300 Hz, in order to suppress their influence on the atan2 detector.
At Fs= ~192 kS/s such a filter can become rather expensive computationally.
Do you think K⁻¹(f) is intentionally unnecessary for the PLL path, or am I missing some mechanism that already makes the LPF inside the loop sufficient?
-
In your proposed structure, K⁻¹(f) is applied only in the S' path for L-R reconstruction, while the PLL path still operates on the original K(f)-shaped signal:
composite → mixdown → LPF → atan2 → PLL
Wouldn't this leave the PLL phase detector exposed to the same altered carrier/sideband ratio introduced by K(f)?
My concern is that without K⁻¹(f) before the PLL path, the close-in sidebands around the subcarrier would still remain relatively elevated with respect to the carrier. In that case it seems that the PLL branch would still require a very narrow LPF, about ~100-300 Hz, in order to suppress their influence on the atan2 detector.
Do you think K⁻¹(f) is intentionally unnecessary for the PLL path, or am I missing some mechanism that already makes the LPF inside the loop sufficient?
The LPF is designed to overcompensate the K(f) shape. Design goals are:
1.Overcompensation of K(f) at any frequency (in order to obtain a signal that does no longer suffer from over-modulation artifacts)
2. Sufficient attenuation for 15 kHz mono components (say -100 dB)
At Fs= ~192 kS/s such a filter can become rather expensive computationally.
Which one do you mean? K-1(z) is a first-oder IIR filter. And LPF can be either a 3rd-order IIR filter, or a FIR block filter (if decimation is used) which is even cheaper than the IIR (only one complex add per sample).
EDIT: See attached plot. With "overcompensated" I mean that the orange curve is completely below the blue one, for any f > 0.
K-1(f) compensates exactly, by reducing far-out sidebands (relative to the subcarrier). And LPF reduces them even more.
[ If LPF had a flat-top shape, then the orange curve could be above the blue one, at low frequencies, before rolling-off steeply. And that's not what we want. ]
-
Thanks, I think I now understand your idea more clearly.
So the LPF after the mixer is not just a simple anti-alias / noise filter, but is actually designed as a combined preconditioning stage that takes the K(f) shaping into account, effectively compensating it in the phase detector path while also ensuring the atan2 input never enters an over-modulation regime.
In other words, the LPF becomes part of the phase detector dynamics (since it follows the VCO-controlled mixer), effectively replacing the need for an explicit K⁻¹(f) stage.
If I understand correctly, this means the loop is intentionally designed so that the LPF effectively counteracts the K(f) shaping in terms of amplitude response, while still preserving sufficient phase information for locking. This leads to a very elegant simplification of the overall demodulation chain, since the K(f) compensation and PLL conditioning are merged into a single stage instead of being handled separately.
One question I still have: since this LPF becomes part of the phase detector transfer function inside the loop, wouldn’t its additional group delay and frequency-dependent phase shift directly affect loop stability and phase margin, especially for very low loop bandwidths (~5–10 Hz)?
-
Thanks, I think I now understand your idea more clearly.
So the LPF after the mixer is not just a simple anti-alias / noise filter, but is actually designed as a combined preconditioning stage that takes the K(f) shaping into account, effectively compensating it in the phase detector path while also ensuring the atan2 input never enters an over-modulation regime.
In other words, the LPF is part of the phase detector dynamics (since it follows the VCO-controlled mixer), so it effectively replaces the need for an explicit K⁻¹(f) stage.
If I understand correctly, this means the loop is intentionally designed so that the LPF dominates the K(f) shaping in terms of amplitude response, while still preserving sufficient phase information for locking.
One question I still have: since this LPF becomes part of the phase detector transfer function inside the loop, wouldn’t its additional group delay and frequency-dependent phase shift directly affect loop stability and phase margin, especially for very low loop bandwidths (~5–10 Hz)?
An additional aim of LPF is the elimination of mono components and negative frequencies in the (down-mixed) composite signal. The input to the phase detector must be strictly limited to the modulated subcarrier. Symmetrical side bands don't matter (as long as they are no over-modulated). As I said, with a loop BW < 10 Hz, the loop's phase margin is stlll acceptable. LPF bandwidth is still much wider than the loop BW. Yet another advantage of the LPF inside the loop is that the VCO will directly track the subcarrier phase, w/o phase delay.
-
This leads to a very elegant simplification of the overall demodulation chain, since the K(f) compensation and PLL conditioning are merged into a single stage instead of being handled separately.
Only for the PLL. You still need to apply K-1(f) to the demoduated S'.
-
[ If LPF had a flat-top shape, then the orange curve could be above the blue one, at low frequencies, before rolling-off steeply. And that's not what we want. ]
Thanks, I think I now understand the issue you are pointing out regarding the LPF shape.
I am currently using a classic windowed Kaiser FIR design for the LPF, but I now see that this approach is essentially "flat-top in behavior" with respect to the K(f) constraint, and therefore not optimal for ensuring the required overcompensation condition.
It makes sense that the filter should not just approximate an ideal low-pass response, but instead be shaped such that it always stays below the inverse K(f) envelope, rather than crossing it at any frequency.
What kind of LPF structure would you recommend in this case?
Would you consider a low-order IIR with explicitly placed poles/zeros, or perhaps a frequency-weighted FIR design to better match this constraint?
-
[ If LPF had a flat-top shape, then the orange curve could be above the blue one, at low frequencies, before rolling-off steeply. And that's not what we want. ]
Thanks, I think I now understand the issue you are pointing out regarding the LPF shape.
I am currently using a classic windowed Kaiser FIR design for the LPF, but I now see that this approach is essentially "flat-top in behavior" with respect to the K(f) constraint, and therefore not optimal for ensuring the required overcompensation condition.
It makes sense that the filter should not just approximate an ideal low-pass response, but instead be shaped such that it always stays below the inverse K(f) envelope, rather than crossing it at any frequency.
What kind of LPF structure would you recommend in this case?
Would you consider a low-order IIR with explicitly placed poles/zeros, or perhaps a frequency-weighted FIR design to better match this constraint?
I already added to a previous message:
% phase detector pre-filter
[b3,a3] = bilinear(1, [tau^3, 3*(tau^2), 3*tau, 1], 1/Fs);
Same tau as K(f).
Or a Hann window as FIR filter (in conjunction with down-sampling).
-
Thanks, this clarifies the idea of using a tau-matched 3rd-order critically damped IIR as a phase detector prefilter.
One question I still have about this approach - how do you evaluate the PLL stability margin in this configuration in practice?
Since the LPF is effectively part of the phase detector transfer function (and also interacts with K(f)-shaping), the open-loop behavior seems to differ significantly from a classical PLL model.
Also, do you take the additional group delay of the 3rd-order filter into account explicitly when tuning the loop bandwidth?
-
Also, do you take the additional group delay of the 3rd-order filter into account explicitly when tuning the loop bandwidth?
At 4.85 Hz gain-crossover frequency (10 Hz loop BW), we talk about ~3ms, or ~5° margin reduction. Don't worry ;)
My experimental PLL locks nicely at 25 Hz offset.
-
Thanks, this clarifies it.
I have another question regarding the phase detector pre-filter concept, but now for CCIR stereo.
Unlike OIRT, CCIR has no K(f) shaping, and there is also a guard region around the 19 kHz pilot, so under ideal conditions there should not be strong close-in sidebands around the carrier itself.
However, in a practical implementation the pilot path is still not perfectly isolated. In my current implementation I use a BPF (and I am considering replacing it with an LPF approach similar to your suggestion). To save computation cost I am not using extremely high attenuation, so some residual high-frequency components from M or S may still leak through, together with noise and other nearby spectral components.
Would a phase detector pre-filter still be beneficial in that case as a way to improve atan2 robustness, or is it mainly useful because of the K(f) shaping in OIRT?
And if it is still beneficial for CCIR, how would you choose its time constant?
Also, I tried adding K(f) shaping to my OIRT test signal and implemented K^-1(f) de-shaping in my OIRT stereo demodulator. Surprisingly, everything works quite well even with my previous PLL implementation: the loop does not lose lock and the lock metric remains stable.
On a test OIRT composite recording with speech the result became noticeably better: bass content reappeared, the stereo effect became stronger, and the sound subjectively became softer and more spacious.
On music recordings, however, some frequency distortions still remain. This may simply be distortion already present in the original recording, although at the moment I cannot completely exclude the possibility that some of it may also come from small PLL jitter caused by strong sideband energy entering the phase detector path.
I will try the new PLL approach and compare the results.
Attached the new synthetic test composite signal with OIRT Polar stereo that now includes K(f) shaping, generated in Octave.
At the moment I am still not completely sure that my implementation is fully correct. The subcarrier level looks reasonably close to the reference OIRT off-air recordings (around -20 dBFS), so that part seems plausible. Although the subcarrier appears to be roughly 6 dB stronger in my file, there is no reliable absolute level reference available because the original WFM signal was not captured for those off-air recordings. The reference recordings were taken directly from the FM discriminator output, therefore only relative amplitude comparisons can be made.
What still puzzles me is that the audio component itself ends up at a rather low level, which is also visible directly in the composite waveform. In real off-air recordings the composite signal appears significantly stronger.
My current suspicion is that perhaps M=(L+R) and S=(L−R) should not actually be divided by 2 in practice. However, the standard explicitly defines them that way, so at the moment I feel I may still be missing something.
-
I have another question regarding the phase detector pre-filter concept, but now for CCIR stereo.
Unlike OIRT, CCIR has no K(f) shaping, and there is also a guard region around the 19 kHz pilot, so under ideal conditions there should not be strong close-in sidebands around the carrier itself.
An atan phase detector cannot deal with out-of-band interfering signals.
So you must ensure sufficient attenuation at f <= 15 kHz and f >= 23 kHz (i.e. outside the guard interval).
For a -100 dB aim, you could try a 3rd order Chebychev2:
[b, a] = cheby2(3, 100, 4000/Fs*2);
Max. group delay is ~1.6 ms, so even shorter than the 3 ms of the 3rd-order critically damped one with tau ~ 1 ms.
[ Btw, the flat top does not matter here -- no need to be strictly below K-1(f). ]
Also, I tried adding K(f) shaping to my OIRT test signal and implemented K^-1(f) de-shaping in my OIRT stereo demodulator. Surprisingly, everything works quite well even with my previous PLL implementation: the loop does not lose lock and the lock metric remains stable.
If the PLL did work without K(f) shaping at the transmitter side, then it will certainly continue to work with transmitter side K(f) shaping. The problem is rather that the receiver side PLL could fail without transmitter side K(f) shaping if the L-R signal happens to be strong and contains very low frequencies (say, a 30 Hz tone at 100% level).
Attached the new synthetic test composite signal with OIRT Polar stereo that now includes K(f) shaping, generated in Octave.
At the moment I am still not completely sure that my implementation is fully correct. The subcarrier level looks reasonably close to the reference OIRT off-air recordings (around -20 dBFS), so that part seems plausible. Although the subcarrier appears to be roughly 6 dB stronger in my file, there is no reliable absolute level reference available because the original WFM signal was not captured for those off-air recordings. The reference recordings were taken directly from the FM discriminator output, therefore only relative amplitude comparisons can be made.
What still puzzles me is that the audio component itself ends up at a rather low level, which is also visible directly in the composite waveform. In real off-air recordings the composite signal appears significantly stronger.
Your composite signal utilizes the [-1, 1] range fully, and sub-carrier level is fine, too. Maybe there are just a few loud peaks and mostly low-volume passages? It also looks like the signal is rather "narrow" (at least on average). Radio stations use dynamic compressors to adjust volume dynamically.
My current suspicion is that perhaps M=(L+R) and S=(L−R) should not actually be divided by 2 in practice. However, the standard explicitly defines them that way, so at the moment I feel I may still be missing something.
If L and R are arbitrary signals with values bounded within the [-1, 1] range (after applying pre-emphasis), then (L+R)/2 and (L-R)/2 and the 0.8 scaling factor make perfectly sense, because they guarantee that the resulting composite signal is strictly bounded within within the [-1, 1] range as well.
-
For a -100 dB aim, you could try a 3rd order Chebychev2:
[b, a] = cheby2(3, 100, 4000/Fs*2);
Max. group delay is ~1.6 ms, so even shorter than the 3 ms of the 3rd-order critically damped one with tau = 1 ms.
[ Btw, the flat top does not matter here -- no need to be strictly below K-1(f). ]
Or, if you want a FIR filter:
b = kaiser(round(0.0010677*Fs), 13); b /= sum(b); a = 1;
Group delay is even lower => 0.54 ms.
-
I think I may have found the reason why the reference off-air recording sounds louder than my synthetic test signal.
The spectrum clearly shows that the off-air recording has a sharp audio cutoff around 10 kHz, while my synthetic test utilize the full 15 kHz audio bandwidth.
Considering the effect of pre-emphasis, limiting the audio bandwidth to 10 kHz allows a noticeably higher input audio level without overdriving the FM modulator. The reason is the rapid increase of composite signal amplitude as the audio content approaches the upper end of the baseband spectrum.
For example, with the pre-emphasis filter, the gain at 10 kHz is approximately 10.362 dB, while at 15 kHz it reaches about 13.656 dB. Therefore, restricting the audio bandwidth to 10 kHz would theoretically allow the input signal level to be increased by about:
13.656 - 10.362 = 3.294 dB
while maintaining the same peak modulation limits.
The attached screenshots illustrate this quite clearly. In the first spectrum view, the audio cutoff can be seen clearly from the sharp boundary of the M component in the real off-air station recording.
The second screenshot also shows how the high-frequency portion of a test sweep (0 - 15 kHz) effectively limits the allowable input signal amplitude. As the sweep approaches the upper end of the audio band, the pre-emphasis increasingly boosts those components, causing them to dominate the composite peak level.
-
I think I may have found the reason why the reference off-air recording sounds louder than my synthetic test signal.
The spectrum clearly shows that the off-air recording has a sharp audio cutoff around 10 kHz, while my synthetic test utilize the full 15 kHz audio bandwidth.
Considering the effect of pre-emphasis, limiting the audio bandwidth to 10 kHz allows a noticeably higher input audio level without overdriving the FM modulator. The reason is the rapid increase of composite signal amplitude as the audio content approaches the upper end of the baseband spectrum.
For example, with the pre-emphasis filter, the gain at 10 kHz is approximately 10.362 dB, while at 15 kHz it reaches about 13.656 dB. Therefore, restricting the audio bandwidth to 10 kHz would theoretically allow the input signal level to be increased by about:
13.656 - 10.362 = 3.294 dB
while maintaining the same peak modulation limits.
Sure, if you statically scale for the worst case pre-emphasis (15 kHz tone with 100% level), then it will become rather quiet on average.
Multiband compressors/limiters even shape the spectrum dynamically on demand.
-
I have revised the 3rd-order IIR pre-filter for the OIR PLL phase detector.
The limiting value for tau3 occurs where the normalized magnitudes |K⁻¹(f)| and |H_LPF(f)| have identical curvature at the peak.
This threshold is tau3 >= sqrt(0.32) * tau, rather than the baseline tau from K(f). [ Thanks to Google for the help with the calculation. ]
Furthermore, attenuation at 16250 Hz is over 106 dB, meaning this condition is also fulfilled.
Max. group delay is only 0.86 ms.
tau3 = sqrt(0.32) * tau;
[b3,a3] = bilinear(1, [tau3^3, 3*(tau3^2), 3*tau3, 1], 1/Fs);
-
I am still fighting with signal level recovery in my FM stereo multiplex decoder for both CCIR and OIRT standards, and I currently have two related issues.
1) OIRT carrier-dependent gain normalization
In OIRT stereo there is an inherent dependency on the 31.25 kHz subcarrier amplitude, because after synchronous AM demodulation of the S signal (L-R component) I need to remove the residual DC component caused by the transmitted carrier.
After thinking about the problem, I came up with the idea of normalizing both M and S channels using the recovered subcarrier level. This allows me to normalize the gain of both components first, and then simply subtract 1.0 as the DC bias from the AM-demodulated S signal.
My reasoning comes from the OIRT modulation formula itself:
composite = 0.8*M + filter(b,a,1+0.8*S) .* real(lo)
where a and b are the K(f) filter coefficients.
If you look carefully at the structure of this equation, it implicitly assumes that after applying the inverse K^-1(f) filter on the receiver side, the recovered carrier amplitude should become exactly 1.0, matching the absolute carrier level originally used at the transmitter. From that perspective, using the recovered subcarrier amplitude as a normalization reference for both M and S initially seemed quite reasonable to me.
Here is my current code for OIRT decoder:
// <...> M = composite sample
// synchronous AM demod S = 2 * real(composite .* conj(lo));
S = 2.0f * M * lo.real();
// inverse K(f)
iirfilt_rrrf_execute(s->inv_kf, S, &S);
// normalize M and S gain by sub-carrier
iirfilt_rrrf_execute(s->lpf_dcS, S-1.0, &G);
G = G + 1.0f; // estimated gain
if (G >= 0.5f && G < 1.5f) {
M = M / G;
S = S / G;
// remove bias DC
S = S - 1.0f;
} else {
S = 0.0f;
}
Here M is the composite signal input. After this stage both M and S are passed through a 15 kHz LPF, then L/R are reconstructed through the stereo matrix, followed by IIR de-emphasis.
At first glance this approach works reasonably well and, importantly, produces noticeably fewer LF transients compared to independently removing DC from the M and S paths.
The only issue is that on real broadcasts the estimated gain G still slightly fluctuating even with a 0.5 Hz cutoff in lpf_dcS. I suspect this happens because I currently use a simple 1st-order LPF with a very gentle roll-off.
One idea is to estimate the carrier level directly inside the PLL instead of after the K^-1(f) equalization filter. However, in that case the recovered carrier amplitude would still include the original K(f) shaping, so I would need to hardcode a compensation factor corresponding to K^-1(0) = +13.979 dB at DC.
I am not sure which approach is considered more correct:
- estimate carrier amplitude directly in the PLL without K^-1(f), which is computationally preferable,
- or estimate it after the inverse K(f) filtering during S recovery.
I also have some doubts whether normalizing both M and S by the recovered subcarrier level is theoretically correct in the first place.
A similar approach could also be applied to CCIR stereo by normalizing M and S using the 19 kHz pilot level. In CCIR this is mostly a gain calibration issue rather than a decoding issue, because DSB-SC has no transmitted carrier and therefore does not require DC removal after demodulation.
My main concern with this approach is that a stereo decoder is expected to continue operating correctly even when the OIRT carrier or CCIR pilot is absent, falling back gracefully to mono operation. With my current normalization scheme, however, the AGC relies entirely on the recovered carrier/pilot amplitude, so in the absence of the subcarrier reference there is nothing left to anchor the gain normalization.
This also implies that if the carrier or pilot disappears, the audio level may change abruptly because the gain normalization mechanism would no longer have a valid reference and would effectively stop operating.
This raises the question of whether normalizing the M and S paths using the recovered carrier/pilot level is actually a good design choice in the first place. If such normalization should be avoided, then what would be the preferred or more correct method of removing the carrier-induced DC bias from the OIRT S channel after demodulation?
2) Correct output gain after de-emphasis
This problem affects both CCIR and OIRT decoders.
The decoder output gain must compensate for the attenuation introduced by the de-emphasis filter. With standard tau=50us, the recovered audio after de-emphasis ends up attenuated by approximately -13.656 dB. This occurs because the FM deviation limit prevents the signal level at the FM modulator input from exceeding unity amplitude. As a result, the transmitter is forced to attenuate the L and R levels to compensate for the gain introduced by the pre-emphasis filter.
However, in practice this attenuation appears to depend on the station. For example, if a station limits its audio bandwidth to around 10 kHz, the effective attenuation after pre-emphasis can be about 4 dB smaller.
As a result, if the receiver blindly applies a fixed +13.656 dB compensation gain, audio from such a 10 kHz-limited station can clip by roughly 4 dB.
So my question is: what is the correct/common way to handle this in a proper FM stereo decoder?
Of course the user can manually adjust gain per station, but that feels like a workaround. Ideally the receiver should recover correct audio levels automatically using only FM deviation and the de-emphasis time constant. How is this normally solved in professional or broadcast-grade receivers?
-
Thanks to the guard bands, the CCIR pilot can be easily isolated, allowing its amplitude to be extracted as an amplitude reference.
In contrast, the OIRT subcarrier cannot be separated from low-frequency sidebands (bass tones) as easily.
If the lpf_dcS filter is not narrow enough, the residual bass tones will unavoidably modulate the AGC gain G(t).
However, in practice this attenuation appears to depend on the station. For example, if a station limits its audio bandwidth to around 10 kHz, the effective attenuation after pre-emphasis can be about 4 dB smaller.
It depends not only on the station but also on the program material itself, meaning a fixed bandwidth limit cannot be expected. Stations heavily utilize dynamic limiters and compressors, often relying on dedicated radio audio processors (Orban Optimod and Telos Omnia seem to be the two market leaders). Their aim is not necessarily an exact mathematical reproduction but rather a subjective, "good-sounding" result, while strictly preventing over-modulation. Ultimately, it is as much an art as it is a science. To deliver a highly dense and loud audio signal that never sounds distorted or harsh to the human ear, these processors rely on a sophisticated combination of psychoacoustic modeling, complex multi-band compressors, and intelligent clippers. While these processors heavily automate the optimization process, they remain highly adjustable, allowing station engineers to sculpt their signature sound.
-
Interestingly, despite the theoretical concerns, the AGC-based approach actually gives noticeably cleaner results on real off-air recordings compared to using separate HPF stages for DC removal in the M and S paths. In particular, it seems to produce significantly fewer LF artifacts and transients.
What still bothers me, however, is how such a scheme should properly handle carrier loss. If the OIRT carrier disappears or becomes unreliable, I would like to avoid audible clicks, sudden gain jumps, or abrupt loudness changes caused by the normalization loop losing its reference.
And what is really frustrating me at the moment is that I still have to adjust the decoder output gain manually, because in practice it appears to depend not only on FM deviation and the de-emphasis time constant, but also quite strongly on the station’s audio bandwidth and processing chain.
As a result, different stations end up requiring different compensation gains, which makes manual per-station adjustment rather inconvenient. At this point I am not even sure what the correct approach should be here, since ideally the decoder should recover proper audio levels automatically without requiring user tuning.
-
It depends not only on the station but also on the program material itself, meaning a fixed bandwidth limit cannot be expected. Stations heavily utilize dynamic limiters and compressors, often relying on dedicated radio audio processors (Orban Optimod and Telos Omnia seem to be the two market leaders). Their aim is not necessarily an exact mathematical reproduction but rather a subjective, "good-sounding" result, while strictly preventing over-modulation. Ultimately, it is as much an art as it is a science. To deliver a highly dense and loud audio signal that never sounds distorted or harsh to the human ear, these processors rely on a sophisticated combination of psychoacoustic modeling, complex multi-band compressors, and intelligent clippers. While these processors heavily automate the optimization process, they remain highly adjustable, allowing station engineers to sculpt their signature sound.
That actually explains a lot and matches what I am observing experimentally. I initially assumed that after compensating for de-emphasis attenuation there should exist some more or less deterministic output gain derived from deviation and the tau constant alone. But in practice modern broadcast chains clearly violate that assumption.
If stations are using aggressive multiband processing, psychoacoustic loudness optimization, adaptive clipping, and dynamically varying bandwidth limits, then the effective gain after pre-emphasis is no longer fixed.
So now I am wondering what would be the best practical approach for a decoder implementation. Would it make more sense to expose the gain correction purely as a manual user option, while by default leaving the signal attenuated after de-emphasis with no compensation at all? Or is it still preferable to apply some automatic correction gain by default — and if so, what would be the most reasonable reference point for choosing it?
-
Different audio materials have distinct crest factors, and human perception of loudness is a matter of acoustic power, and also frequency-dependent. Scaling different audio sources to the same peak voltage does not result in equal perceived loudness. To implement a perceptual audio AGC, you should instead normalize the K-weighted average power rather than the voltage. [ Or the K-weighted average power sum, for multi-channel audio. ]
Seems there also exists an ITU Recommendation
Algorithms to measure audio programme loudness and true-peak audio level (https://www.itu.int/dms_pubrec/itu-r/rec/bs/R-REC-BS.1770-5-202311-I!!PDF-E.pdf)
While radio processors/compressors can (and certainly do) reduce the crest factor, full equalization would likely compromise perceived sound quality too much. Given the strict constraint of maximum FM deviation (which bounds the composite signal's peak level) it is IMO practically impossible to transmit different audio sources such that a perfectly uniform perceived loudness is achieved after demodulation.
I don't think that classical analog broadcast receivers implemented audio AGC, let alone a perceptual one. Their AGC was limited to the RF/IF stages to stabilize signal reception. On the other hand, in an SDR, implementing this feature has a low cost-overhead, as it is handled entirely in software.
EDIT: However, it should be noted that real-time audio AGC unavoidably introduces additional dynamic range compression, even if it operates with a longer time constant than typical compressors. If a station intentionally reduces compression to retain a higher dynamic range (e.g. for classical music), a receiver-side AGC would compromise that effort.
-
Interestingly, despite the theoretical concerns, the AGC-based approach actually gives noticeably cleaner results on real off-air recordings compared to using separate HPF stages for DC removal in the M and S paths. In particular, it seems to produce significantly fewer LF artifacts and transients.
The question is, what's the actual defect of the composite signal? Do you refer to the two records you uploded to a file share?
-
Yes, I am referring to those two signals.
I believe the issue is simply that their amplitude is not normalized, since they were recorded directly from the discriminator output of a receiver without any reference to an absolute level. As a result, the “1.0” amplitude in these files does not correspond to the maximum FM deviation at the discriminator input.
Because of this, without carrier-based AGC the reconstructed M and S levels are incorrectly scaled. Once I apply normalization based on the recovered subcarrier, the result becomes significantly more consistent.
So my current assumption is that the problem is mainly due to missing absolute level calibration in the recorded signal rather than a flaw in the demodulation itself.
However, the idea of normalizing the M and S amplitudes using the recovered carrier still seems relevant, as it provides a convenient way to quickly restore consistent gain for both components. But it is not entirely clear whether this is a robust approach for a stereo decoder in general, because there may be hidden pitfalls that make it unsuitable in some cases.
Nevertheless, I find this method appealing because it allows simultaneous gain correction for both M and S, preventing them from drifting apart in level, and it also effectively removes the DC offset without requiring separate per-channel filtering stages.
In essence, I replaced two separate filtering stages for M and S with a single unified gain normalization stage, which also implicitly removes the DC offset introduced by DSB demodulation.
Previously, both M and S components were drifting around zero and were stabilized using independent LPF structures. Now this instability is essentially gone, and both signals behave as stably as in a CCIR decoder.
The only remaining challenge is reliably isolating the carrier level from the adjacent sideband components, which is non-trivial in the OIRT case. At the moment I am experimenting with deriving the carrier amplitude directly from a PLL implemented with IIR filtering. In that case, the carrier level could be extracted directly from the phase detector after frequency translation to baseband, which might provide a cleaner and more convenient reference.
PS: For the past couple of days I was trying to understand the source of a spurious crosstalk between channels in my CCIR decoder, at a level around -90 dB. Today I finally found the cause.
It turned out to be related to a recently introduced mechanism for smooth decoder engagement based on PLL lock status, where the S signal was being attenuated under conditions of weak PLL lock. Because the lock indicator was not strictly binary, this effectively introduced partial gain modulation of the S path, which in turn resulted in unintended channel leakage.
After fixing this behavior, the issue disappeared, and the CCIR stereo separation is now reaching around 140 dB. :)
-
Yes, I am referring to those two signals.
I believe the issue is simply that their amplitude is not normalized, since they were recorded directly from the discriminator output of a receiver without any reference to an absolute level. As a result, the “1.0” amplitude in these files does not correspond to the maximum FM deviation at the discriminator input.
Because of this, without carrier-based AGC the reconstructed M and S levels are incorrectly scaled. Once I apply normalization based on the recovered subcarrier, the result becomes significantly more consistent.
Normalized amplitude is not required. Pure scaling of the composite signal affects M and S equally, so it does not impact stereo decoding (except for the absolute loudness of the resulting L and R channels). However, these two recordings have defects that go beyond mere overall scaling. For example, a significant and unexpected phase variation/modulation of the subcarrier is clearly evident. I don't know the root cause. Maybe it is distortion from the IF filter or discriminator? Maybe multipath effects? Or perhaps just noise? :-//
EDIT: Do you know what kind of radio was used for the two recordings? Old analog Soviet radio receiver?
-
However, these two recordings have defects that go beyond mere overall scaling. For example, a significant and unexpected phase variation/modulation of the subcarrier is clearly evident. I don't know the root cause. Maybe it is distortion from the IF filter or discriminator? Maybe multipath effects? Or perhaps just noise? :-//
EDIT: Do you know what kind of radio was used for the two recordings? Old analog Soviet radio receiver?
Yes. As @Njk mentioned in the thread (https://www.eevblog.com/forum/rf-microwave/wfm-stereo-broadcasting-standards-oirt-vs-ccir/msg6035139/#msg6035139), the recordings were made using a Soviet analog radio receiver, the "Leningrad RP-015 (https://www.radiomuseum.org/r/lenin_nova_leningrad_015_c.html)". He also posted a fragment of the schematic there, showing that the recording was taken from test point KT7 ~80 мВ (located near transistor V32, see slightly above and to the left of the center of the schematic image).
(https://www.eevblog.com/forum/rf-microwave/wfm-stereo-broadcasting-standards-oirt-vs-ccir/?action=dlattach;attach=2656979;image)
So the signal was not captured from a calibrated SDR or measurement setup, but rather from an internal test point in the receiver's FM stereo demodulation chain. That makes me wonder whether some of the anomalies we are seeing could indeed be related to the receiver itself rather than the original broadcast signal. For example, imperfections in the IF filter, discriminator, limiter stages, or other parts of the analog signal path could potentially introduce the phase and amplitude artifacts observed in the recordings.
-
Of course, having a proper IQ recording made with an RTL-SDR or similar receiver would make the analysis much easier.
Unfortunately, OIRT stereo broadcasts are quite rare nowadays. They can mostly still be found in regions where the 65.9–74 MHz FM band remains in use, but not every station in that band transmits stereo, and not all stereo stations necessarily use the OIRT stereo standard.
I do have an IQ recording of a Belarusian FM 70.43 MHz station with OIRT stereo received in France during a DX propagation event, but the recording quality is far from ideal. The signal is unstable and periodically fades into the noise, making it difficult to draw reliable conclusions from it. However, the recording is still good enough to estimate the level of the 31.25 kHz subcarrier.
If anyone lives in an area where OIRT FM broadcasting is still active and has access to a receiver capable of making IQ recordings, I would greatly appreciate a recording of an OIRT stereo station. Such a recording would be extremely valuable for validating and improving the OIRT polar stereo DSP decoder.
-
I extracted a short segment from an IQ recording of a Belarusian OIRT FM 70.43 MHz station "Kultura" with polar stereo, received in France during a DX propagation event. An interesting detail is that this FM signal was received on 70.43 MHz from a transmitter located about 2000 km away. :)
Below is a small Octave script for loading the file, FM demod, displaying MPX composite spectrum, extract and playing the mono L+R component as an example.
pkg load signal;
[wfm,fs] = audioread('test-oirt-70430kHz-sample.flac');
if size(wfm, 2) ~= 2
error('IQ stream expected');
end
wfm = wfm(:,1) + 1j*wfm(:,2);
% FM demod using polar discriminator
mpx = angle(wfm(2:end) .* conj(wfm(1:end-1)));
mpx = mpx * fs / (2*pi*50000); % normalize for FM deviation 50 kHz
fprintf('FM demod done\n');
% LPF 96 kHz/-100 dB & decimate to fs=192 kHz
[n, Wn, beta, ftype] = kaiserord([87000, 96000], [1, 0], [ 0.0057564, 10^(-100/20) ], fs);
n += mod(n,2); % force even order
lpf96k = fir1(n, Wn, ftype, kaiser(n+1,beta));
mpx = filter(lpf96k,1, mpx);
fprintf('LPF 96k done\n');
D = round(fs/192000);
mpx = downsample(mpx, D);
fs = fs / D;
%audiowrite('test-oirt-mpx.wav', mpx, fs, 'BitsPerSample', 24);
if 1
% plot MPX composite spectrum
NFFT = min(2^20, 2^(nextpow2(length(mpx))-1));
f = (0:NFFT/2) * fs / NFFT;
w = kaiser(NFFT, 10.0613); % sidelobe -100 dB
h = fft(mpx(1:NFFT) .* w);
h = h(1:NFFT/2+1);
mag = abs(h) / sum(w);
mag = 20*log10(max(mag, eps));
figure;
plot(f/1000, mag);
xlim([min(f) max(f)]/1000);
ylim([-90 10]);
grid on; grid minor; set (gca, 'gridalpha', 0.2); set(gca, 'minorgridalpha', 0.05);
xlabel('kHz');
ylabel('dB');
title(sprintf('FFT %g kpts, fs=%g kHz', NFFT/1024, fs/1000));
end
if 1
% de-emphasis 50us
tau = 50e-6;
wcw = 2.0*fs*tan((1.0/tau)/(2.0*fs)); % prewarped wc=1/tau
z1 = -1.0; % zero
p1 = (2.0*fs - wcw)/(2.0*fs + wcw); % pole
k = (1.0 - p1) / (1.0 - z1); % gain
b = [ 1, -z1 ] * k;
a = [ 1, -p1 ];
M = filter(b,a, mpx);
% LPF 15 kHz/-80 dB & decimate to fs=48 kHz
[n, Wn, beta, ftype] = kaiserord([15000, 31250-15000], [1, 0], [ 0.0057564, 10^(-80/20) ], fs);
n += mod(n,2); % force even order
lpf15k = fir1(n, Wn, ftype, kaiser(n+1,beta));
M = filter(lpf15k,1, M);
fprintf('LPF 15k done\n');
D = round(fs / 48000);
M = downsample(M, D);
fs = fs / D;
sound(M, fs, 16);
end
Despite the bad SNR quality, the advantage of this record sample is that it contains the original FM signal in its raw form, without any intermediate processing. This makes it possible to directly observe the actual level of the 31.25 kHz subcarrier in a real-world OIRT broadcast signal.
Another advantage of this sample is that it contains a real off-air FM signal with genuine periodic fading caused by changing DX propagation conditions. As a result, it is quite useful for testing and tuning PLL behavior under realistic reception conditions.
-
To be honest, it’s a pretty horrible signal. ;) The peak of the largest fading hump at ~6.5s only gets around 10 dB SNR. The other humps are even weaker, and the SNR goes completely negative in the fading valleys. Since the SNR peaks are right on the edge of (or even below) the FM mono threshold, I'm not surprised that the PLL can just barely track the subcarrier on the hills, but has no chance during the deep fades.
EDIT: Attached plot for estimated carrier/noise ratio.
-
Overall, I tested the triple-RC IIR filter you suggested, and it appears to work well for PLL.
tau = 6.4 / (2*pi*1000);
tau3 = sqrt(0.32) * tau;
%[b,a] = bilinear(1, [tau3^3, 3*(tau3^2), 3*tau3, 1], 1/fs);
% H(s) = 1 / (tau3*s + 1)^3;
% H(z) = [1 + 3*z^-1 + 3*z^-2 + z^-3] / [A^3 + 3*A^2*B*z^-1 + 3*A*B^2*z^-2 + B^3*z^-3];
K = tau3*2.0*fs;
A = 1.0+K;
B = 1.0-K;
A2 = A*A;
A3 = A2*A;
B2 = B*B;
B3 = B2*B;
b = [ 1.0/A3, 3.0/A3, 3.0/A3, 1.0/A3 ];
a = [ 1.0, 3.0*B/A, 3.0*B2/A2, B3/A3 ];
I have also integrated it into the PLL loop of my stereo decoder. Based on my tests, the preliminary BPF at the PLL input has almost no noticeable effect on performance, so I removed it entirely and switched the PLL input to a real-valued signal to reduce processing overhead.
So far, everything seems to be working well — both OIRT and CCIR stereo signals decode properly. With this IIR PLL lock is much more stable. For CCIR, I am using the same IIR filter as well.
An additional benefit is that the real-time performance of the CCIR decoder has improved significantly, with processing speed nearly doubling. Now on raspi 4 it eats just 3% CPU load in realtime for 240 kHz IQ bandwidth, include WFM demod and resampling 240->192->48 :)
By the way, could you please explain the purpose of the sqrt(0.32) factor used in the equation of tau3?
As I understand it, this IIR is effectively a cascade of three identical RC low-pass sections, and the coefficient is intended to compensate for the bandwidth reduction compared to a single RC filter. If so, was that the rationale behind introducing it?
What is not entirely clear to me is whether this is a standard design approach for a three-stage RC filter. Where does this coefficient come from, and how is it derived?
-
EDIT: Attached plot for estimated carrier/noise ratio.
here is how this file works with my PLL with your tripple-RC IIR filter in the loop:
-
By the way, could you please explain the purpose of the sqrt(0.32) factor used in the equation of tau3?
As I understand it, this IIR is effectively a cascade of three identical RC low-pass sections, and the coefficient is intended to compensate for the bandwidth reduction compared to a single RC filter. If so, was that the rationale behind introducing it?
What is not entirely clear to me is whether this is a standard design approach for a three-stage RC filter. Where does this coefficient come from, and how is it derived?
This is a third-order critically damped IIR filter *). A filter of at least third-order is necessary to achieve sufficient attenuation for the mono signal below 15 kHz. Setting τ₃ = τ · √0.32 yields an equal curvature of │H(f)│ and │K⁻¹(f)│ at their peak, ensuring that H(f) overcompensates for K(f). Credits go to Google's math skills for identifying the ratio of √0.32, but my verification of the frequency response proves that it successfully delivers the intended overcompensation.
EDIT: *) Critically damped, to get a round top and prevent a flat top response.
EDIT: τ₃ = τ · √0.32 is the minimum value. You can narrow the filter further as long as the PLL tolerates the increased group delay—a narrower filter will also yield a lower ENBW.
-
For CCIR, I am using the same IIR filter as well.
Since you only have a 4 kHz guard band on each side, the filter for the 19 kHz pilot PLL should be narrower, to prevent leakage from f < 15 kHz and from f > 23 kHz.
A tau3 = 0.0018464 should result in -100 dB at +-4 kHz.
I'd still use a critically damped one, as it also has a lower ENBW than filters with a flat top.
EDIT: However, this is also a matter of PLL phase margin. A ~5.5 ms group delay should still work with a 5 Hz loop bandwidth and damping 0.707. With a 10 Hz loop bandwidth, phase margin becomes borderline. A Chebyshev Type II filter like cheby2(3, 100, 4000/(Fs/2)) would reduce the group delay to only ~2.3 ms, but it allows more noise to pass through. [ To clarify, since PLL "loop bandwidth" lacks a universal definition, I am specifically referring to the closed-loop single-sided noise bandwidth. ]
-
Here is a small status update on the OIRT stereo decoder.
1) I finally identified the factor that was limiting channel separation. The issue turned out to be insufficient suppression of the DC component in the S signal. To estimate the DC level, I replaced the single EMA LPF with a cascaded structure consisting of three EMA low-pass filters with fc = 5 Hz.
This change completely resolved the poor stereo separation issue and improved channel separation from approximately 60 dB to 110 dB. Now OIRT separation is about ~110 dB and CCIR separation is > 140 dB :) Unfortunately OIRT decoder is 3 times slower than CCIR, since it requires to apply K^-1(f) de-shaping IIR filter and DC estimation to decode S component.
I experimented with several filter orders and configurations, but consistently good results were obtained only with a cascade of three EMA stages. The frequency response of the currently used filter for DC estimation is shown below.
This made me consider whether it might be more appropriate to implement the IIR LPF used in the PLL as a cascade of three EMA sections with an equivalent time constant tau3?
2) During debugging SOS filter for DC estimation, another question arose regarding accurate group delay estimation for cascaded sequence of IIR filters. Conventional numerical approaches tend to become unstable for such filters, especially in highly attenuated stopband regions. Currently I am using the following approach to estimate it:
fs = 192000; % sample rate
nsos = 3; % number of second-order sections (sos)
b = [...
0.000163611228344962, 0, 0; ...
0.000163611228344962, 0, 0; ...
0.000163611228344962, 0, 0; ...
];
a = [...
1, -0.99983638525009155, 0; ...
1, -0.99983638525009155, 0; ...
1, -0.99983638525009155, 0; ...
];
for i=1:nsos
if abs(sum(b(i,:))) < 1e-6 fprintf('warn: SOS %d: sum(b)==%g\n', i, sum(b(i,:))); end
if abs(sum(a(i,:))) < 1e-6 fprintf('warn: SOS %d: sum(a)==%g\n', i, sum(a(i,:))); end
end
nfft=2^nextpow2(fs);
w = 2*pi*(0:nfft-1)/nfft;
z = exp(1j*w);
z1 = z.^-1;
z2 = z.^-2;
z3 = z.^-3;
H = ones(1,nfft);
gd = zeros(1, nfft);
for i = 1:nsos
b0 = b(i,1); b1 = b(i,2); b2 = b(i,3);
a0 = a(i,1); a1 = a(i,2); a2 = a(i,3);
% H(z)
B = b0 + b1*z1 + b2*z2;
A = a0 + a1*z1 + a2*z2;
H = H .* (B ./ A);
% derivatives
dB = -(b1*z2 + 2*b2*z3);
dA = -(a1*z2 + 2*a2*z3);
% group delay
gd = gd + real( (z .* dA ./ A) - (z .* dB ./ B) );
end
f = w/(2*pi) - 0.5;
H = fftshift(H);
gd = fftshift(gd);
H = 20*log10(max(abs(H),eps));
f = (f * fs)/1000; xlabel_txt='kHz';
fprintf('DC: %.3f dB\n', H(find(f==0)));
figure;
subplot(2,1,1); hold on;
plot(f,H, '-','Color',[0 0 1], 'LineWidth',2);
axis([0 0.3 -120 10]);
grid on; grid minor; set (gca, 'gridalpha', 0.2); set(gca, 'minorgridalpha', 0.05);
xlabel(xlabel_txt);
ylabel('dB');
subplot(2,1,2); hold on;
plot(f,gd, '-','Color',[0 0 1], 'LineWidth',2);
axis([0 0.3 -1000 +20000]);
grid on; grid minor; set (gca, 'gridalpha', 0.2); set(gca, 'minorgridalpha', 0.05);
xlabel(xlabel_txt);
ylabel('Group delay [samples]');
Does this look like a correct method for estimating group delay in the sequence of several second-order IIR sections?
3) I also discovered an interesting detail related to FM modulation.
After fixing a gain error in the (M,S) → (L,R) reconstruction matrix, I noticed that real FM broadcast stations appear to feed the FM modulator directly from the output of the pre-emphasis filter without any subsequent level normalization. If this assumption is correct, the actual frequency deviation exceeds the nominal 50/75 kHz value by approximately the gain of the pre-emphasis filter, which is roughly 14–17 dB depending on the selected time constant tau.
At first I was skeptical about this conclusion, so I checked several local CCIR FM stations. In all cases, the pilot tone amplitude was around 0.09 for an assumed deviation of 75 kHz. This appears reasonable, since the standard pilot injection level is 9–10%.
However, with this deviation, the amplitude of the M signal (L+R), after applying the 1/0.9 correction factor (since the M component occupies 90% of the total deviation according to the standard), exceeds the nominal full-scale level by approximately 17 dB. Interestingly, this value is very close to the gain introduced by the pre-emphasis filter itself.
This led me to suspect that many real-world transmitters may not normalize the signal level after pre-emphasis, effectively operating with peak deviations significantly above the nominal 50/75 kHz value due to the pre-emphasis boost.
My question is: how well does this align with actual broadcasting standards and industry practice? Is this something specific to the local stations I measured, or is it generally accepted practice in FM broadcasting systems?
-
Hmmm, something isn't plausible here. At a 384 kSa/s sampling rate, the angle(cross product) FM demodulator is not able to "see" 17 dB above 75 kHz peak deviation, because its ±180° output range limits the detectable deviation to ±192 kHz (which is only 8 dB above 75 kHz) before a phase slip/wrap-around happens. Furthermore, a peak deviation of ±530.96 kHz would far exceed regulatory limits, and I also cannot imagine that hardware receivers possess such a wide IF bandwidth.
EDIT: Can you post (say) 10s of the I/Q IF signal?
-
Maximum group delay tg_max of a first order lowpass = 1 / (2 * pi * fc), and it occurs at f = 0.
If you chain 3 of them, you get tg_max = 3 / (2 * pi * fc).
For arbitrary filters, I usually calculate the group delay numerically: I use [H,f] = freqz(b, a, fs/2, fs), which results in a 1 Hz frequency spacing to obtain the complex frequency response H. Then, I compute the group delay by taking the negative numerical derivative of the unwrapped phase: gd = -diff(unwrap(angle(H))) / (2 * pi * delta_f). At a 1 Hz frequency spacing, delta_f = 1.
In Octave, there is also a grpdelay() function, but at least in my version, it gives strange results.
-
Maximum group delay tg_max of a first order lowpass = 1 / (2 * pi * fc), and it occurs at f = 0.
If you chain 3 of them, you get tg_max = 3 / (2 * pi * fc).
Thanks, I think your formula refers to the analog filter H(s) = 1 / (1+tau*s)^3. Is that correct?
For the digital EMA IIR filter H(z) = alpha / (1-beta*z^-1), where beta = exp(-1/(fs*tau)), alpha = 1 - beta, the DC group delay should be:
gd(0) = beta / alpha
For fc=5 Hz and fs=192 kHz, this gives a DC group delay of about 6111 samples for a single EMA stage, or about 18333 samples for three cascaded stages... I hope I estimate it the right way...
For arbitrary filters, I usually calculate the group delay numerically: I use [H,f] = freqz(b, a, fs/2, fs), which results in a 1 Hz frequency spacing to obtain the complex frequency response H. Then, I compute the group delay by taking the negative numerical derivative of the unwrapped phase: gd = -diff(unwrap(angle(H))) / (2 * pi * delta_f). At a 1 Hz frequency spacing, delta_f = 1.
Yes, the numerical method
gd = -diff(unwrap(angle(H))) / (2*pi*delta_f);
does work, and in fact that was my starting point. However, I ran into several issues:
- one sample is lost because of diff();
- very large spikes appear near deep response minima due to numerical sensitivity;
- results near the frequency range boundaries become unreliable and usually show random huge spikes.
I found that
w = 2*pi*(0:nfft-1)/nfft;
gd = -gradient(unwrap(angle(H)), w);
behaves somewhat a little bit better because it preserves the number of points, but it still suffers from boundary artifacts and numerical instability in regions where the magnitude response becomes very small.
Eventually I switched to the analytical derivative approach similar to the one used in liquid-dsp. As I understand, it is based on the group-delay identity:
\$\tau_{g}\left( \omega \right)=-\Re \left(
z\frac{H'(z)}{H(z)}
\right)\$
which avoids phase unwrapping and numerical differentiation entirely, and can be evaluated directly from the numerator and denominator polynomials. In my tests, this method is much more stable near deep stopband nulls and at the frequency range boundaries. However, this approach is new to me, so I'm not yet completely sure that I am applying it correctly.
If I understand correctly, this identity follows from the definition of group delay:
\$\tau_{g}\left( \omega \right)=-\frac{d\phi(\omega)}{d\omega}\$
by expressing the phase derivative through the logarithmic derivative of the transfer function.
One of the main difficulties when analyzing IIR filters, especially cascaded or higher-order, is the large numerical dynamic range involved. The frequency response can easily span hundreds of dB, and intermediate values may approach the limits of double-precision floating-point arithmetic. This is particularly noticeable when evaluating group delay in deep stopband regions, where the magnitude response becomes extremely small.
That is the reason why I started looking for more numerically robust analysis methods. My goal is not to improve the filter itself, but rather to obtain more reliable estimates of its characteristics (magnitude, phase and group delay) when conventional numerical differentiation begins to break down.
In Octave, there is also a grpdelay() function, but at least in my version, it gives strange results.
Regarding grpdelay(), it actually seems to work reasonably well and is more robust than the diff() or gradient() approaches. My main issue is that it returns results only on the positive-frequency interval 0..+pi, while I prefer to analyze the full spectrum -pi...+pi, especially when comparing different filter structures.
-
Thanks, I think your formula refers to the analog filter H(s) = 1 / (1+tau*s)^3. Is that correct?
For the digital EMA IIR filter H(z) = alpha / (1-beta*z^-1), where beta = exp(-1/(fs*tau)), alpha = 1 - beta, the DC group delay should be:
gd(0) = beta / alpha
For fc=5 Hz and fs=192 kHz, this gives a DC group delay of about 6111 samples for a single EMA stage, or about 18333 samples for three cascaded stages... I hope I estimate it the right way...
Well, not really much difference to the analog one:
>> 18333/192000
ans = 0.095484
>> 3/(2*pi*5)
ans = 0.095493
That is the reason why I started looking for more numerically robust analysis methods. My goal is not to improve the filter itself, but rather to obtain more reliable estimates of its characteristics (magnitude, phase and group delay) when conventional numerical differentiation begins to break down.
I didn't dive deeply into analytic calculation yet. Occational large spikes are frequently not real, but just a numerical unwrap problem which can be solved by adding or subtracting 2*pi from the affected diff(angle(H))(index). And at zeros, the angle() is of course undefined, as atan2(0,0) is ambiguous and could be any angle. The function just returns zero by definition. So don't expect a reasonable numerical gd at/near zeros.
Regarding grpdelay(), it actually seems to work reasonably well and is more robust than the diff() or gradient() approaches. My main issue is that it returns results only on the positive-frequency interval 0..+pi, while I prefer to analyze the full spectrum -pi...+pi, especially when comparing different filter structures.
For filters with real coefficients, H(-f) is simply a conjugate complex mirror of H(f). The double-sided frequency response becomes only relevant if the filter has complex coefficients. Do you really need it?
EDIT: Seems that grpdelay() accepts a 'whole' parameter to calculate the group delay between 0 and 2*pi.
But as said, mine seems to be broken anyway (older version).
-
I didn't dive deeply into analytic calculation yet. Occational large spikes are frequently not real, but just a numerical unwrap problem which can be solved by adding or subtracting 2*pi from the affected diff(angle(H))(index). And at zeros, the angle() is of course undefined, as atan2(0,0) is ambiguous and could be any angle. The function just returns zero by definition. So don't expect a reasonable numerical gd at/near zeros.
Yes, I agree regarding atan2(0,0). At an exact zero of the transfer function the phase is undefined, so we cannot expect a meaningful numerical group delay there. But what I find more problematic in practice is that with real filters we are usually not dealing with exact zeros, but with extremely small values that are merely close to zero. In those regions numerical precision starts to dominate the result. Small errors in the complex response can produce large phase fluctuations, and after differentiation these fluctuations may turn into significant noise or huge spikes in the estimated group delay.
I have also noticed that for some cascaded IIR filters, even though the filter itself is perfectly usable in practice, functions such as freqz() or grpdelay() may occasionally return Inf or otherwise suspicious values at isolated frequencies. This seems to be more related to numerical conditioning than to the actual filter behavior.
What I like about the polynomial-derivative approach is that it appears to be much less sensitive to these numerical issues. At least in my experiments it produces smooth and consistent group-delay curves even in regions where the conventional phase-differentiation approach becomes noisy or unstable. Of course I am still evaluating it and trying to understand its limitations.
For filters with real coefficients, H(-f) is simply a conjugate complex mirror of H(f). The double-sided frequency response becomes only relevant if the filter has complex coefficients. Do you really need it?
EDIT: Seems that grpdelay() accepts a 'whole' parameter to calculate the group delay between 0 and 2*pi.
But as said, mine seems to be broken anyway (older version).
Thanks for mentioning the 'whole' option in grpdelay(). I was not aware of that parameter. I tried it on a very simple SOS structure consisting of three cascaded EMA IIR filters. Unfortunately, in my case grpdelay() produces a large number of warnings such as: "warning: grpdelay: setting group delay to 0 at singularity"
As a result, the computed group-delay curve is essentially pinned to zero and does not show the expected large delay near DC.
In contrast, the analytical derivative approach remains well-behaved and produces a realistic group-delay estimate near DC.
The same issue happens with freqz() when I use 4 EMA sections in series - the response obtained with freqz() near DC looks like garbage.
Since the filter itself is simply a cascade of three stable first-order EMA sections, I think this is a numerical conditioning issue with freqz and grpdelay() rather than an actual singularity of the filter.
Here is octave test code:
pkg load signal;
fs = 192000; % sample rate
fc = 5; % cut-off frequency
nsos = 3; % number of second-order sections (sos)
nfft = 2^nextpow2(fs);
% EMA IIR section
beta = exp(-2*pi*fc/fs);
alpha = 1 - beta;
b = [alpha, 0, 0];
a = [1, -beta, 0];
% build SOS (nsos cascaded identical sections)
b = repmat(b, nsos, 1);
a = repmat(a, nsos, 1);
% sanity check
for i=1:nsos
if abs(sum(b(i,:))) < 1e-6 fprintf('warn: SOS %d: sum(b)==%g\n', i, sum(b(i,:))); end
if abs(sum(a(i,:))) < 1e-6 fprintf('warn: SOS %d: sum(a)==%g\n', i, sum(a(i,:))); end
end
%-----------------------------------------------------------------------
% Polynomial derivative method
w = 2*pi*(0:nfft-1)/nfft;
z = exp(1j*w);
z1 = z.^-1;
z2 = z.^-2;
z3 = z.^-3;
H1 = ones(1,nfft);
gd1 = zeros(1, nfft);
for i = 1:nsos
b0 = b(i,1); b1 = b(i,2); b2 = b(i,3);
a0 = a(i,1); a1 = a(i,2); a2 = a(i,3);
% H(z)
B = b0 + b1*z1 + b2*z2;
A = a0 + a1*z1 + a2*z2;
H1 = H1 .* (B ./ A);
% derivatives
dB = -(b1*z2 + 2*b2*z3);
dA = -(a1*z2 + 2*a2*z3);
% group delay
gd1 = gd1 + real( (z .* dA ./ A) - (z .* dB ./ B) );
end
% Convert to FFT-style axis: -fs/2 ... +fs/2
H1 = fftshift(H1);
gd1 = fftshift(gd1);
f1 = w/(2*pi) - 0.5; % normalized frequency -0.5..+0.5
%-----------------------------------------------------------------------
% freqz / grpdelay method
B = 1;
A = 1;
for i = 1:nsos
B = conv(B, b(i,:));
A = conv(A, a(i,:));
end
[H2,w] = freqz(B, A, nfft, 'whole');
[gd2,wg] = grpdelay(B, A, nfft, 'whole');
if max(abs(w-wg)) > 0
warning('frequency vectors are not identical!');
end
% Convert to FFT-style axis: -fs/2 ... +fs/2
H2 = fftshift(H2);
gd2 = fftshift(gd2);
f2 = w/(2*pi) - 0.5; % normalized frequency -0.5..+0.5
%-----------------------------------------------------------------------
% Plot
% unit scale
f1 = f1 * fs / 1000; % kHz
f2 = f2 * fs / 1000; % kHz
H1 = 20*log10(max(abs(H1), eps)); % dB
H2 = 20*log10(max(abs(H2), eps)); % dB
figure;
subplot(2,1,1); hold on;
plot(f1, H1, 'Color',[0 0 1], 'LineWidth',3, 'DisplayName','Polynomial H(z) eval');
plot(f2, H2, 'Color',[1 0 0], 'LineWidth',1, 'DisplayName','freqz()');
axis([-0.1 0.1 -80 10]);
grid on; grid minor; set (gca, 'gridalpha', 0.2); set(gca, 'minorgridalpha', 0.05);
xlabel('kHz');
ylabel('dB');
legend();
title(sprintf('%d cascaded second-order sections', nsos));
subplot(2,1,2); hold on;
plot(f1, gd1, 'Color',[0 0 1], 'LineWidth',3, 'DisplayName','Polynomial derivative');
plot(f2, gd2, 'Color',[1 0 0], 'LineWidth',1, 'DisplayName','grpdelay()');
axis([-0.1 0.1 -1000 +20000]);
grid on; grid minor; set (gca, 'gridalpha', 0.2); set(gca, 'minorgridalpha', 0.05);
xlabel('kHz');
ylabel('Group delay [samples]');
legend();
just try to change nsos from 1 to 4 and you will see the difference.
-
Unfortunately, in my case grpdelay() produces a large number of warnings such as: "warning: grpdelay: setting group delay to 0 at singularity"
As a result, the computed group-delay curve is essentially pinned to zero and does not show the expected large delay near DC.
In contrast, the analytical derivative approach remains well-behaved and produces a realistic group-delay estimate near DC.
The "singularity" issue is exactly my problem with the grpdelay() function.
OTOH, numerical derivatives worked fine in my case, and analytic ones are certainly expected to work fine as well.
EDIT: If freqz is inaccurate with direct form, you can also apply it to each sos stage separately and multiply the resulting H(f).
For a low cut-off frequency like 5 Hz, a 1 Hz step is likely too wide for a reasonable approximation; you rather need 0.1 Hz or even less.
-
Hmmm, something isn't plausible here. At a 384 kSa/s sampling rate, the angle(cross product) FM demodulator is not able to "see" 17 dB above 75 kHz peak deviation, because its ±180° output range limits the detectable deviation to ±192 kHz (which is only 8 dB above 75 kHz) before a phase slip/wrap-around happens. Furthermore, a peak deviation of ±530.96 kHz would far exceed regulatory limits, and I also cannot imagine that hardware receivers possess such a wide IF bandwidth.
EDIT: Can you post (say) 10s of the I/Q IF signal?
Your reasoning is generally sound, and I agree that ±530.96 kHz looks strange. However, I still observe a clear mismatch in amplitude scaling in my measurements, which does not fully align with the expected deviation limits.
This suggests that there might be some hidden factors involved. I will capture a clean IQ recording of a strong FM broadcast station and share it so we can analyze it more precisely in a controlled way.
-
One thing that can be mistaken for overmodulation is noise.
If you have a (say) 20 dB RF SNR with Gaussian noise, the demodulator turns that into frequency deviation noise with a standard deviation of roughly 6.1 kHz. Six-sigma spikes look exactly like ~37 kHz overmodulation spikes. And at 384 kSa/s, six-sigma events are expected every 22 minutes on average.
And with only 10 dB SNR, standard deviation is about 20 kHz, and six-sima spikes are about 120 kHz.
Still it cannot explain 17 dB.
-
I'm currently looking for a suitable FM station to record, preferably one with high SNR and no significant adjacent-channel activity or interference within ±200 kHz, so that the measurements are not affected by interference.
Interestingly, it appears that I catch a station while its FM deviation was changing in real-time. The effect is quite noticeable even by ear: when the deviation increases to nearly twice its normal value, the DSP chain begins to exhibit obvious overload artifacts. As you can see on the spectrum, peak deviation is about ±200 kHz in this extended mode, so the total station bandwidth is more than 400 kHz! :o And this station doesn't have RDS.
What is particularly interesting is that this increased deviation seems to occur specifically during commercial breaks. Based on several observations, the station appears to switch to a significantly higher modulation level whenever advertisements are broadcast.
The screenshot shows the moment when the station switches back from the increased-deviation mode to its normal deviation level.
-
Here is a IF recording sample of the same station operating with its normal (lower-deviation) modulation mode.
https://gofile.io/d/PKoeKF
added my measurements for this file. M amplitude is shown without de-emphasis (just 15 kHz LPF). The pilot tone amplitude at the beginning of the graph smoothly reaches a level, because it is passed through an LPF with a cutoff of 1 Hz to cut off the AC components around the pilot.
Most of these peaks above 1 are high-frequency components amplified by the pre-emphasis filter. I think this confirms that the pre-emphasis filter gain is not compensated on the transmitter side.
It appears that the magic here is how FM modulation responds to excessive amplitude in the high-frequency components of the modulating signal.
-
I have run into another issue.
The CCIR/OIRT stereo encoder/decoder utility includes an optional FM modulator/demodulator stage, allowing an FM-modulated IF stream to be decoded directly. During testing, I noticed that when the encoder output is fed directly into the decoder, enabling the FM mod => FM demod chain introduces faint but clearly audible whistle-like tones in the recovered audio. If the FM modulation/demodulation stage is bypassed entirely, the output audio is clean.
Initial testing showed that increasing the sample rate used by the FM modulator from 384 kHz to 768 kHz does not improve the situation. For FM modulation and demodulation I am currently using the liquid-dsp freqmod_* and freqdem_* functions. Interestingly, I observe essentially the same behavior when performing FM modulation in Octave using cumtrapz, which probably suggests that the issue is not specific to the liquid-dsp implementation.
At this point the root cause remains unclear. However, the tests performed so far indicate that the M (L+R) component remains essentially unaffected, while the artifacts appear in the S (L−R) component. This suggests that the distortion mechanism is somehow interacting with the stereo-difference path rather than the mono component.
Do you have any ideas about possible sources of the noise and distortion observed in the S (L−R) signal due to adding FM mod => FM demod chain?
It is worth noting that this issue is not specific to OIRT stereo decoding. The same behavior is present with CCIR stereo as well, where recovery of the S component is comparatively straightforward and consists essentially of frequency shift of the DSB-SC subchannel to baseband followed by low-pass filtering.
-
If you modulate using
phase = cumsum(75000 * mpx * 2 * pi * Ts);
fm = exp(1i * [0; phase]);
then
mpx_demod = angle(fm(2:end) .* conj(fm(1:end-1))) / 75000 / 2 / pi / Ts;
reconstructs mpx exactly (withn the limits of numerical precision).
>> max(abs(mpx_demod - mpx))
ans = 1.1657e-14
Note that neither cumsum() nor cumtrapz() are perfect integrators. However, cumsum() is the inverse operation of diff(), therefore the angle difference demodulator reconstructs mpx exactly if the FM modulator uses cumsum() for the integration.
EDIT: Compared to cumsum(), cumtrapz() attenuates high frequencies slightly, as if a [0.5, 0.5] FIR kernel were applied.
Note that H_cumtrapz(z) = (0.5 + 0.5z⁻¹) / (1 - z⁻¹) can be decomposed into (0.5 + 0.5z⁻¹) * H_cumsum(z).
The first term (0.5 + 0.5z⁻¹) is a 2-tap FIR filter, and H_cumsum(z) = 1 / (1 - z⁻¹).
The initial conditions (i.e., the initial carrier phase) differ as well, but this is not really relevant for FM.
Also note that an FM signal has theoretically infinite bandwidth. If you cut-off at say +-100 kHz, then the reconstruction is no longer exact.
-
Here is a IF recording sample of the same station operating with its normal (lower-deviation) modulation mode.
https://gofile.io/d/PKoeKF (https://gofile.io/d/PKoeKF)
The ETSI mask (https://www.etsi.org/deliver/etsi_en/302000_302099/302018/02.01.01_60/en_302018v020101p.pdf) requires the maximum out-of-band power density to drop from 0 dBc/kHz at a 100 kHz offset to -80 dBc/kHz at a 200 kHz offset, measured with a 1 kHz reference bandwidth and max hold. I simulated measurement conditions with an STFT using a 770 point Blackman-Harris window (ENBW ~1 kHz). The SDR capture is quite noisy, with the max hold noise floor sitting around -40 dBc/kHz, making it hard to assess directly. However, with an apparent slope of only roughly -0.4 dB/kHz, I find it unlikely that the signal will achieve the required -80 dBc/kHz attenuation at the 200 kHz offset.
When demodulating the signal, noise can explain part of the observed overmodulation, but definitively not all of it.
I also think the signal is unlikely to be compliant.
-
The ETSI mask (https://www.etsi.org/deliver/etsi_en/302000_302099/302018/02.01.01_60/en_302018v020101p.pdf) requires the maximum out-of-band power density to drop from 0 dBc/kHz at a 100 kHz offset to -80 dBc/kHz at a 200 kHz offset, measured with a 1 kHz reference bandwidth and max hold.
Thanks for the reference to the ETSI spectral mask for FM broadcasting — that’s useful context.
However, the requirement that the mask already starts effectively constraining power at a 100 kHz offset looks a bit inconsistent with the nominal FM deviation of 75 kHz. Here is why I find it questionable.
In a standard CCIR stereo multiplex signal, the composite baseband is limited by the upper edge of the DSB-modulated S component, i.e. 38 kHz + 15 kHz = 53 kHz. So the highest significant baseband frequency is about 53 kHz.
Applying Carson’s rule, the occupied FM bandwidth can be approximated as:
BW ≈ 2 (Δf + f_max)
With Δf = 75 kHz and f_max ≈ 53 kHz, this gives:
BW ≈ 2 × (75 kHz + 53 kHz) ≈ 256 kHz, i.e. roughly ±128 kHz around the carrier.
This seems to directly conflict with a ETSI spectral mask that already requires strong suppression at more than ±100 kHz offsets. 🤔
I wonder how this can be explained? After all, the spectral mask essentially only allows for ±100 kHz carrier deviation, while the CCIR stereo composite requires ±128 kHz. 🤔
I simulated measurement conditions with an STFT using a 770 point Blackman-Harris window (ENBW ~1 kHz). The SDR capture is quite noisy, with the max hold noise floor sitting around -40 dBc/kHz, making it hard to assess directly. However, with an apparent slope of only roughly -0.4 dB/kHz, I find it unlikely that the signal will achieve the required -80 dBc/kHz attenuation at the 200 kHz offset.
When demodulating the signal, noise can explain part of the observed overmodulation, but definitively not all of it.
I also think the signal is unlikely to be compliant.
Yes, I also noticed that the SNR is only around ~40 dB, which is not particularly high for this kind of spectral measurement. Still, this is actually one of the stronger stations in terms of received power in my setup. The only station with noticeably higher level is about 10 dB stronger, but it is heavily contaminated with spurs, so it is not really suitable for analysis.
What is interesting is that most stations in my environment are even lower, often around ~30 dB SNR, yet they still decode perfectly well and behave like normal FM broadcast signals in a typical urban reception scenario.
Another interesting observation is that this particular station shows a relatively moderate occupied bandwidth compared to others. Visually, it appears to fit reasonably well within the ETSI spectral mask. In contrast, several other stations exhibit noticeably wider spectra, which is likely not only due to RDS, but probably also due to differences in audio processing and modulation chain behavior.
One more curious point is that this station has been observed to temporarily expand its occupied bandwidth by nearly a factor of two, reaching approximately ±200 kHz during live broadcast. This happens specifically during commercial breaks. It strongly suggests that the advertising insertion path is locally generated and not processed through the same audio chain as the main program.
My current hypothesis is that this bandwidth expansion is caused by the absence of the usual broadcast processing chain (in particular, dynamic compression and multi-band limiting) during local ad insertion, which leads to significantly higher peak levels and occasional overmodulation.
I’ve come up with a few ideas for testing the FM signal bandwidth and need to experiment with them. I’ll add the results once the experiments are completed.
-
Yes, I also noticed that the SNR is only around ~40 dB, which is not particularly high for this kind of spectral measurement.
M2 M4 SNR Estimation gives only ~21 dB SNR for this IF signal.
[ The roughly -40 I mentioned was not the SNR, but the relative power density (PSD) at the noise shoulder (per kHz BW), measured using max hold over ~2 ms frames with a 1 kHz ENBW. ]