1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
| close all; clear all; clc;
% ****************************
% (1) CREATE A SIGNAL + NOISE
% ****************************
% Generate some binary data for our "signal"
Nbits = 8192;
data = randi([0,1],Nbits,1);
% Let's modulate the data as QPSK and say the baud rate is 10kHz.
syms = bi2de(reshape(data,2,[]).');
x = exp(1j*2*pi*syms/4);
fbaud = 1e4;
% Upsample to 10x baud rate
upsamp = 10;
y = resample(x,upsamp,1); % (let Matlab use its default interpolation filter)
fsamp_y = fbaud * upsamp;
% Add some noise. Let's give it about 1% of the signal power (SNR ~20dB).
Psig = y'*y / length(y);
n = sqrt(0.01 * Psig) * complex(randn(size(y)), randn(size(y))) / sqrt(2);
z = y + n;
fsamp_z = fsamp_y;
% Print the actual SNR in dB.
SNRdb = 10*log10((y'*y) / (n'*n))
% Plot the QPSK constellation
figure();
plot(z(1:upsamp:end),'x');
grid on; axis equal;
xlabel('Real'); ylabel('Imag')
title('Noisy QPSK constellation');
% Compute spectrum
Nfft = 1024;
Pz = pwelch(z,hann(Nfft),[],[],fsamp_z);
% Unwrap spectrum
Pz = fftshift(Pz);
f = [ -(ceil((Nfft-1)/2):-1:1), 0, (1:floor((Nfft-1)/2)) ] * fsamp_z / Nfft;
% Plot spectrum
figure();
plot(f, 10*log10(Pz));
axis tight; grid on;
xlabel('Freq (Hz)'); ylabel('PSD (dB)');
title('Spectrum of 10x oversampled signal');
% *************************
% (2) DECIMATE AND COMPARE
% *************************
% Decimate and retain whatever default filter, h, is used
downsamp = 4;
[zdec, h] = resample(z,1,downsamp);
fsamp_zdec = fsamp_z / downsamp;
% Compute the frequency response of the decimating filter
H = freqz(h,1,Nfft,"whole",fsamp_z);
% Unwrap spectrum
H = fftshift(H);
% Compute spectrum of decimated signal with the same resolution as before
Nfft_zdec = Nfft/downsamp;
Pzdec = pwelch(zdec,hann(Nfft_zdec),[],[],fsamp_zdec);
% Unwrap spectrum
Pzdec = fftshift(Pzdec);
fzdec = [ -(ceil((Nfft_zdec-1)/2):-1:1), 0, (1:floor((Nfft_zdec-1)/2)) ] * fsamp_zdec / Nfft_zdec;
figure(); hold on;
ax = plotyy(f, 10*log10(Pz), f,10*log10(abs(H)));
plot(fzdec, 10*log10(Pzdec), '--g');
axis tight;
grid on;
legend('Decimated Signal + Noise', 'Input Signal + Noise', 'Filter');
line([-1,-1]*fsamp_zdec/2,ylim,'LineStyle','--','Color','k');
line([1,1]*fsamp_zdec/2,ylim,'LineStyle','--','Color','k');
xlabel('Freq (Hz)'); ylabel(ax(1), 'PSD (dB)'); ylabel(ax(2), '(dB)'); |