Low-pass filter in carrier modulation and demodulation

Question

I'm designing a project in which an array is passed through quadrature amplitude modulation (QAM) modulator, and then do carrier modulation, make it playable with the sound() command, then demodulate it back for QAM demodulation.

Firstly, I have used the standard way of QAM modulation:

M = 16;
x = randint(5000, 1, M);
y = modulate(modem.qammod(M), x);

Then, I wrote my own carrier modulation function:

function [out] = carriermodulation(x)
fs = 16000;
T = 1.0 / 4000;
fc = 8000;
Q = real(x);
I = imag(x);
t = 0:T:(size(x))*T;
C1 = zeros(size(x), 1);
C2 = zeros(size(x), 1);
for i = 1:size(x)
    C1(i) = I(i)*sin(2*pi*(fc)*t(i));
    C2(i) = Q(i)*sin(2*pi*fc*t(i) + pi/2);
end

out = C1 + C2;

No problem so far. But when I was done with my demodulation function, I found that the result is different from the original value (the QAM modulator output).

function [out] = carrierdemodulation(x)
fs = 16000;
T = 1.0 / 4000;
fc = 8000;
t = 0:T:(size(x))*T;
A1 = zeros( size(x), 1);
A2 = zeros( size(x), 1);
for i = 1:size(x)
    A1(i) = x(i)*sin( 2*pi*(fc)*t(i));
    A2(i) = x(i)*cos( 2*pi*(fc)*t(i));
end
A1 = sqrt(A1);
A2 = sqrt(A2);
out = A1 + A2;

I think my modulation part is right. The only problem I think I have is I don't have a low-pass filter (LPF) for the demodulation. And I should not calculate A1 and A2 directly. How do I to add an LPF to my demodulation code such that the output is the same as the original?

Answer 1

You need a low-pass filter at the receiver after coherent demodulation, that's right. But there's also a problem with your modulation. In your example, the symbol rate Rs is less than the angular carrier frequency w_c which potentially causes overlapping of spectra at the receiver. Consequently, reconstruction of the information signal will be impossible. Additionaly, in your example fc * T = 2. That means that the argument of the sine function is an integer multiple of 2pi and therefore always zero.

What you need is an impulse shaper (can be implemented as a low-pass filter) at the transmitter with bandwidth w_g >= R/2. It should be a so-called Nyquist lowpass. The carrier frequency must satisfy w_c > w_g.

I've written a MATLAB script that does impulse shaping, modulation, demodulation, filtering and sampling, so that the transmitted signal can be reconstructed.

First we define the parameters, create random bits and do the mapping as you've already done. Than a very simple impulse response for impulse shaping is used, namely a rectangular impulse. In the real world we're going from digital to analog domain here, but as this is a computer model, we represent the analog signal by a discrete one with sampling frequency f_s. The impulse shaper is simple, because it just repeats each sample L times.

M = 16; % QAM order
fs = 16000; % Sampling frequency in Hz
Ts = 1/fs; % Sampling interval in s
fc = 1000; % Carrier frequency in Hz (must be < fs/2 and > fg)
Rs = 100; % Symbol rate
Ns = 20; % Number of symbols

x = randint(Ns, 1, M);
y = modulate(modem.qammod(M), x);

L = fs / Rs; % Oversampling factor

% Impulse shaping
y_a = reshape(repmat(y', L, 1), 1, length(y)*L);

Now modulation. I used a carrier frequency that satisfies the above conditions: it's higher than the signal bandwidth and can still be represented with the sampling frequency used.

%% Modulation
I = real(y_a);
Q = imag(y_a);
t = 0 : Ts : (length(y_a) - 1) * Ts;
C1 = I .* sin(2*pi * fc * t);
C2 = Q .* cos(2*pi * fc * t);
s = C1 + C2;

Demodlation is straightforward...

%% Demodulation
r_I = s .* sin(2*pi * fc * t);
r_Q = s .* -cos(2*pi * fc * t);

To remove the spectral tributaries at 2f_c after demodulation a low-pass filter is required. I used the MATLAB FDATool to create the filter and part of the following code. Remember: the signal bandwidth is Rs/2, and the unwanted tributaries begin at 2*fc - Rs/2. This is how Fpass and Fstop are found. (It might be useful to relax these requirements a little bit.)

%% Filter

% Design filter with least-squares method
N     = 50;           % Order
Fpass = Rs/2;         % Passband Frequency
Fstop = 2*fc - Rs/2;  % Stopband Frequency
Wpass = 1;            % Passband Weight
Wstop = 1;            % Stopband Weight

% Calculate the coefficients using the FIRLS function.
b  = firls(N, [0 Fpass Fstop fs/2]/(fs/2), [1 1 0 0], [Wpass Wstop]);

% Filtering
w_I = filter(b, 1, r_I);
w_Q = filter(b, 1, r_Q);

After filtering we still have to sample the received signal. Here it's just a downsampling. I used a phase offset of L/2 to avoid the filter transitions.

%% Sampling
u_I = downsample(w_I, L, L/2);
u_Q = downsample(w_Q, L, L/2);

Finally, plot the constellation diagram and get a nice 16-QAM constellation:

plot(u_I, u_Q, '.');

You can find the complete code here.

Your question touched a lot of topics of both DSP and MATLAB programming. I could not go into much detail everywhere. If you have specific questions about 16-QAM modulation and demodulation the place to be is probably the Stack Exchange site Signal Processing.