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.
A1=x.*sin(2*pi*fc*t);
will result in the sameA1
as your whole loop does. – Dan