Matlab code for SSB demodulation using phasing method

Started by danielot 5 years ago11 replieslatest reply 5 years ago1055 views


I'm trying to simulate the phasing method to perform SSB AM demodulation in Matlab but something is going wrong. I'm using as reference the literature online below and I think I'm following it closely, unless there's a flaw in my understanding (this is probably happening).

My idea is to simply quadrature downconvert my IF (called IF1 in the code) coming from the ADC samples to a lower IF (IF2 in the code) and then apply -45 deg phase shift in the I branch and +45 deg in the Q branch with bandpass filters with identical magnitude response. To get USB I'm then doing I - Q and to get LSB I'm doing I + Q.

I would appreciate if you could take a look in my code below. I would expect that the signals 'usb' and 'lsb' in the end of the code would have one of the two sidebands (signal and image) suppressed, but they are still there after my processing.



% == Signal generation ==
cic_nsecs = 2;  % Number of sections of prototype CIC decimation filter

nperiods = 10 + cic_nsecs;

Mif1 = 98;
N = 383;

t = ((0:N*nperiods-1) - cic_nsecs*N)';

usb = cos(2*pi*(Mif1+0.2)/N*t);        % Upper sideband
lsb = 0.3*cos(2*pi*(Mif1-0.4)/N*t);    % Lower sideband

x = usb + lsb;

% == Build bandpass filters to replace Hilbert transform ==
Mif2 = 8;

% Build prototype CIC decimation filter
lpf_decim = 1;
for i=1:cic_nsecs
    lpf_decim = conv(lpf_decim,ones(1,N)/N);

% Translate LPF CIC decimator to BPF centered at IF2
cos_minus_45 = 2*cos(2*pi*Mif2/N*(0:N*cic_nsecs-1)-pi/4);
cos_plus_45 = 2*cos(2*pi*Mif2/N*(0:N*cic_nsecs-1)+pi/4);

coeff_fir_I = cos_minus_45.*[0 lpf_decim];
coeff_fir_Q = cos_plus_45.*[0 lpf_decim];

% == SSB demodulation phasing method ==
Mnco = Mif1 - Mif2;
nco = exp(1j*(2*pi*Mnco/N*t));

x_if2 = x.*nco;
x_if2_filtered = filter(coeff_fir_I, 1, real(x_if2)) + 1j*filter(coeff_fir_Q, 1, imag(x_if2));

usb = real(x_if2_filtered) - imag(x_if2_filtered);
lsb = real(x_if2_filtered) + imag(x_if2_filtered);

% == Plots ==
y = [x x_if2 x_if2_filtered usb lsb];
y = y(cic_nsecs*N+1:end,:);
npts = size(y,1);

Y = fft(y)/npts;
for i=1:size(Y,2)
    Y(:,i) = fftshift(Y(:,i));

leg = {'x', 'x_{shifted}', 'x_{analytic}', 'x_{only USB}', 'x_{only LSB}'};
f = (-ceil(npts/2):floor(npts/2)-1)/npts*N;

ax(1) = subplot(211);
plot(f, real(Y));
ylabel('real part')
grid on
hold on
ax(2) = subplot(212);
plot(f, imag(Y));
ylabel('imaginary part')
grid on
[ - ]
Reply by Rick LyonsMay 12, 2017

Hello Daniel.  I just saw your post this morning. Have you solved your 'MATLAB SSB demodulation' problem?

[ - ]
Reply by danielotMay 12, 2017

Hi Rick, thanks for your interest.

Yes, I did. I'm pasting my finished "product" in the end of this message (sorry for the formatting, it seems DSPrelated doesn't allow code formatting in the replies).

Instead of designing a Hilbert filter or doing -45 deg / +45 deg phase shift in I/Q branches I stepped back and used the Matlab hilbert() function.

What I've learned from this interaction was that there's a need to have an abrupt 90 deg phase shift in the frequency between both sidebands, requiring a "brickwall cutoff" (I like kaz's exprresion :), which will always make it difficult to recover the carrier (or DC in baseband) with a low latency filter. This is the same tradeoff as in Weaver's method. Do you confirm that?

My application is very sensitive to DC (actually from ~10 uHz to 1 kHz) since it directly translates to synchrotron light source beam stability, and it is also very sensitive to digital processing latency since the system runs in feedback. So, maybe the SSB demod isn't suited for me.

My finished code below. Thanks.

% == Test signals ==
freqstep = 0.1;

nperiods = 1/freqstep;

M = 4;
N = 20;

t = (0:N*nperiods-1)';

x_usb = zeros(size(t));
x_lsb = zeros(size(t));
bb_usb = zeros(size(t));
bb_lsb = zeros(size(t));

for i=1:10
    omega_sidebands = 2*pi*freqstep*i/N*t + pi/4;
    usb_amplitude_profile = (0.1*i + 0.2);
    lsb_amplitude_profile = 0.5*(1.1 - 0.1*i);
    x_usb = x_usb + usb_amplitude_profile*cos(2*pi*M/N*t + omega_sidebands);
    x_lsb = x_lsb + lsb_amplitude_profile*cos(2*pi*M/N*t - omega_sidebands);
    bb_usb = bb_usb + usb_amplitude_profile*cos(omega_sidebands);
    bb_lsb = bb_lsb + lsb_amplitude_profile*cos(-omega_sidebands);

x = [x_usb  x_lsb];
bb = [bb_usb  bb_lsb];

% == SSB demodulation using phasing method ==
nco = exp(1j*(-2*pi*M/N*t));

I = x.*repmat(real(nco),1,2);
Q = x.*repmat(imag(nco),1,2);
Q_hilbert = imag(hilbert(Q));

usb_recovered = I - Q_hilbert;
lsb_recovered = I + Q_hilbert;

% == Plots ==
y = {bb, x, I, Q, Q_hilbert, usb_recovered, lsb_recovered};

for i=1:length(y)
    npts = size(y{i},1);
    Y{i} = fft(y{i})/npts;
    for j=1:size(Y{i},2)
        Y{i}(:,j) = fftshift(Y{i}(:,j));
    f{i} = (-ceil(npts/2):floor(npts/2)-1)/npts*N;

graph_names = {'baseband', 'RF', 'I', 'Q', 'Q_{Hilbert}', 'Recovered USB', 'Recovered LSB'};
leg = {'USB','LSB'};

ax = zeros(2*length(y),1);
j = 1;
close all
for i=1:length(y)
    ax(j) = subplot(211);
    plot(f{i}, real(Y{i}));
    ylabel('real part');
    grid on;
    ax(j+1) = subplot(212);
    plot(f{i}, imag(Y{i}));
    ylabel('imaginary part');
    grid on;
    j = j+2;


[ - ]
Reply by kazMay 11, 2017

Just a thought, my understanding of USB & LSB modulation is that you have to modulate by removing one sideband before transmission. You are not doing that, you just start from x = sum of two frequencies.

[ - ]
Reply by danielotMay 11, 2017

I'm trying to do demodulation, not modulation.

In the end of the processing I would like to recover x with only one of the sidebands. Once I have this I can downconvert to baseband.

Is this possible? I couldn't find any straight forward example of such processing, so maybe I'm understanding SSB demodulation wrongly.

For the SSB modulation, as you mentioned, I found this visual interpreation which is incredibly nice!


[ - ]
Reply by kazMay 11, 2017

I understand you are doing demodulation but you still need to start from a properly SSB modulated signal for your demod testing.

[ - ]
Reply by danielotMay 11, 2017

Well, so maybe this means it is impossible to achieve what I wanted in first place.

The problem I'm trying to solve is the following: I have a carrier at 500 MHz and AM modulation with ~580 kHz bandwidth (form 500 MHz - 290 kHz to 500 MHz + 290 kHz). I want to translate to baseband only one of the sidebands around 500 MHz, for instance, 500 MHz to 500 MHz + 290 kHz, and discard the other one completely.

This is intended for beam position measurement in a particle accelerator and there are some specific reasons why I want to avoid USB to interfere with LSB when downconverting.

[ - ]
Reply by kazMay 11, 2017

Have you thought of downconverting to baseband followed by complex filter?

[ - ]
Reply by danielotMay 11, 2017

Actually not. What kind of complex filtering would you propose? I understand once I complex downconvert to baseband I'll get one sideband from -290 kHz to DC and the other sideband from DC to 290 kHz. If you propose I complex filter out frequencies below DC, I fear I will loose much of the power at DC, which is critical in my application.

I was also thinking of Weaver's method, but in this case too my filter (simple low pass) would need very steep cut off, which will also hurt latency, another critical parameter in my application (the measurement coming out this processing will be used in a fast feedback).

My original code wants to integrate decimation plus Hilbert transform in a single step so as to minimize latency. This would be beautiful provided it works, which is not the case :)

If you have references for the complex filtering you mentioned please let me know. Thanks!

[ - ]
Reply by kazMay 11, 2017

As any filter, yes you will have to cut sharp enough and tolerate latency. You will need really good filter with brickwall cutoff to do what you target.  

You might have to give it more thoughts...

Matlab has complex filter function. Or you can design your own by setting negative bins to zero then getting coeffs.

[ - ]
Reply by kazMay 11, 2017

if you want avoid complex filter you can downconvert somewhere near dc apply ordinary filter on both Re/Im to remove one sideband then downconvert to dc

[ - ]
Reply by danielotMay 11, 2017

Nice, got it!

This last tip seems to bring to Weaver's method, which now I can see is as troublesome as the phasing method when it comes to achieving low latency. A trade-off will have to be made...

Thanks, kaz.