## Matlab code for SSB demodulation using phasing method

Started by 7 years ago11 replieslatest reply 7 years ago1313 views

Hi,

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.

Thanks,

Daniel

% == 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);
end

% 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));
end

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

figure;
ax(1) = subplot(211);
plot(f, real(Y));
ylabel('real part')
grid on
hold on
legend(leg);
ax(2) = subplot(212);
plot(f, imag(Y));
ylabel('imaginary part')
grid on
legend(leg)
[ - ]

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

[ - ]

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);
end

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));
end
f{i} = (-ceil(npts/2):floor(npts/2)-1)/npts*N;
end

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)
figure;
ax(j) = subplot(211);
plot(f{i}, real(Y{i}));
ylabel('real part');
grid on;
legend(leg);
title(graph_names{i});
ax(j+1) = subplot(212);
plot(f{i}, imag(Y{i}));
ylabel('imaginary part');
grid on;
j = j+2;
end

[ - ]

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.

[ - ]

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! http://k6jca.blogspot.com.br/2017/02/sdr-notes-wea...

Thanks!

[ - ]

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

[ - ]

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.

[ - ]

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

[ - ]

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!

[ - ]

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.

[ - ]