实验

实验内容

产生门函数信号

实验代码
clc,clear;
% 参数设定
tao = 2;
fs = 10;
ts = 1/fs;
T = 2;
t = -T:ts:T-ts;
K1 = abs(-T + tao/2)/ts;
K2 = abs(T-tao/2)/ts;
M = tao/ts;
% 产生门函数信号
x = [zeros(1, K1) 2*ones(1, M) zeros(1, K2)];
% 绘图
plot(t, x, 'o-');grid
title('门函数信号波形');
xlabel('t(s)');ylabel('幅度');
实验效果

高斯白噪声通过带通和低通滤波器

实验代码
clc,clear all;
fs = 2000;
ts = 1/fs;
T = 50;
t = 0:ts:T-ts;
N = T*fs;
df = fs/N;
f = [0:df:(N-1)*df]-fs/2;
mean = 0;
n0 = 0.002;
delta2 = n0*fs;
noisein = mean+sqrt(delta2)*randn(1, N);
bin = 1000;
[num1, x1] = hist(noisein, bin);
p1 = num1/N;
dx = (max(noisein)-min(noisein))/bin;
figure(1);
subplot(2, 1, 1);
plot(x1, p1/dx);
xlabel('样值大小');ylabel('取值概率');title('高斯白噪声概率密度');

Xnisein = fftshift(fft(noisein, N)/fs);
Pnoisein = (abs(Xnisein)).^2;
subplot(2, 1, 2);
plot(f,Pnoisein);
xlabel('f/Hz');ylabel('幅度');title('白噪声功率密度');

figure(2);
subplot(3, 1, 1);
bp = [500 600]/(fs/2);
b = fir1(1000, bp);
a = 1;
HP = freqz(b, a, f, fs);
plot(f, abs(HP));
xlabel('f/Hz'),ylabel('幅度');title('带通滤波器传输特性');

subplot(3, 1, 2);
n_bp = filter(b, a, noisein);
[num2, x2] = hist(n_bp(1000:N), bin);
plot(x2, num2/(dx*N));
xlabel('样值大小');ylabel('取值概率');title('带通滤波器输出的瞬时概率密度');

subplot(3, 1, 3);
Xn_bp = fftshift(fft(n_bp, N)/fs);
Pn_bp = (abs(Xn_bp)).^2;
plot(f, Pn_bp);
xlabel('f/Hz'),ylabel('幅度');title('带通滤波器输出的噪声功率谱密度');

figure(3);
subplot(3, 1, 1);
bl = 500/(fs/2);
b = fir1(1000, bl);
a = 1;
HL = freqz(b, a, f, fs);
plot(f, abs(HL));
xlabel('f/Hz'),ylabel('幅度');title('低通滤波器传输特性');

subplot(3, 1, 2);
n_bl = filter(b, a, noisein);
% n_bl = conv(b, noisein)*ts;
[num3, x3] = hist(n_bl(1000:N), bin);
dx = (max(n_bl(1000:N)) - min(n_bl(1000:N)))/bin;
plot(x3, num3/(dx*N));
xlabel('f/Hz'),ylabel('幅度');title('低通滤波器输出的噪声功率谱密度');

subplot(3, 1, 3);
Xn_bl = fftshift(fft(n_bl, N)/fs);
Pn_bl = (abs(Xn_bl)).^2;
plot(f, Pn_bl);
xlabel('f/Hz'),ylabel('幅度');title('低通滤波器输出的噪声功率谱密度');
实验效果

单/双极性不归零与归零码

实验代码
clc,clear;
Tb = 1;
M = 200;
ts = Tb/M;
K = 0.5;
N = 100;
s = randi(1, N);
KK = floor(M*K);

for i = 1:length(s)
    s1((i-1)*M+1:i*M) = s(i)*ones(1, M);
end

for i = 1:length(s)
    s2((i-1)*M+1:i*M) = s(i)*[ones(1,KK) zeros(1,M-KK)];
end

for i = 1:length(s)
    s3((i-1)*M+1:i*M) = (-1)^(s(i)-1)*ones(1, M);
end

for i = 1:length(s)
    s4((i-1)*M+1:i*M) = (-1)^(s(i)-1)*[ones(1,KK) zeros(1,M-KK)];
end

index = 10;
t1 = 0:Tb:index*Tb-Tb;
t2 = 0:ts:index*Tb-ts;
subplot(5, 1, 1);
stem(t1,s(1:index));grid
axis([0 index*Tb -1.2 1.2]);
xlabel('t');ylabel('幅度');title('数字信息');

subplot(5, 1, 2);
plot(t2,s1(1:index*M));grid
axis([0 index*Tb -1.2 1.2]);
xlabel('t');ylabel('幅度');title('单极性不归零码');

subplot(5, 1, 3);
plot(t2,s2(1:index*M));grid
axis([0 index*Tb -1.2 1.2]);
xlabel('t');ylabel('幅度');title('单极性归零码');

subplot(5, 1, 4);
plot(t2,s3(1:index*M));grid
axis([0 index*Tb -1.2 1.2]);
xlabel('t');ylabel('幅度');title('双极性不归零码');

subplot(5, 1, 5);
plot(t2,s4(1:index*M));grid
axis([0 index*Tb -1.2 1.2]);
xlabel('t');ylabel('幅度');title('双极性归零码');
实验效果

正弦函数的均均与量化

实验代码
clc,clear;
T = 1.2;
t = 0:0.01:T;
x = sin(2*pi*t);

fs = 10;
ts = 1/fs;
t1 = 0:ts:T;
xs = sin(2*pi*t1);

n = 4;
x_width = max(x)-min(x);
delta = x_width/n;
xx = min(x):delta:max(x);
q = min(x) + delta/2:delta:max(x);
for i = 1:length(xs);
    if xs(i) < xx(2)
        xq(i) = q(1);
    end
    k = 2;
    while k < n+1;
        if(xs(i) > xx(k)) & (xs(i) <= xx(k + 1))
            xq(i) = q(k);
            k = n+1;
        end
        k = k+1;
    end
end

plot(t, x, 'r');hold on
plot(t1, xs, 'go');hold on
plot(t1, xq, '*');
axis([0 1.2 -1.2 1.2])
title('模拟信源与抽样信号');ylabel('幅度');xlabel('t/sec');legend('模拟信源', '抽样值', '量化信号');
实验效果

窄带高斯随机过程功率谱及包络分布

实验代码
clc,clear;
f0 = 20000;
fs = 100000;
ts = 1/fs;
df = 0.1;
N = fs/df;
n0 = 0.0001;
delta2 = n0*fs;
x1 = sqrt(delta2)*randn(1,N);
x2 = sqrt(delta2)*randn(1,N);
a = [1 -0.9];
b = 1;
xc = filter(b, a, x1);
xs = filter(b, a, x2);
t = 0:ts:(N-1)*ts;
band_pass_process = xc.*cos(2*pi*f0*t)-xs.*sin(2*pi*f0*t);
env = sqrt(xc.^2 + xs.^2);
Xbp = fftshift(fft(band_pass_process, N)/fs);
f = [0:df:(N-1)*df] - fs/2;
plot(f, (abs(Xbp)).^2);
xlabel('f/Hz');ylabel('幅度');title('窄带高斯随机过程功率谱')

figure;
bin = 100;
[sum, x] = hist(env, bin);
dx = (max(env)-min(env))/bin;
penv = sum/(N*dx);
plot(x,penv);
xlabel('x/v');ylabel('取值概率');title('窄带高斯随机过程包络瞬时值的概率密度函数')
实验效果

差分编码仿真及其波形

实验代码
clc,clear;
% 参数设置
Tb = 1;
M = 200;
ts = Tb/M;
K = 0.5;
KK = floor(M*K);
a = [0 1 0 1 1 1 0 0 1 0 1];

% 传号差分码
% 单极性码运算
b(1) = 0;
for i = 1:length(a)
    b(i+1) = xor(a(i), b(i));
end

% 形成指定的波形
for i = 1:length(b)
    b_s_NRZ((i-1)*M+1:i*M) = b(i)*ones(1,M);      % 单极性不归零波形
    b_s_RZ((i-1)*M+1:i*M) = b(i)*[ones(1,KK) zeros(1,M-KK)];   % 单极性归零波形
    b_d_NRZ((i-1)*M+1:i*M) = (-1)^(b(i)-1)*ones(1,M);          % 双极性不归零波形
    b_d_RZ((i-1)*M+1:i*M) = (-1)^(b(i)-1)*[ones(1,KK) zeros(1,M-KK)];   %双极性归零波形
end

% 双极性码运算
aa = (-1).^(a-1);
bb(1) = -1;
for i = 1:length(aa)
    bb(i+1) = -aa(i)*bb(i);
end

% 形成指定的波形
for i = 1:length(bb)
    bb_NRZ((i-1)*M+1:i*M) = b(i)*ones(1,M);
    bb_RZ((i-1)*M+1:i*M) = b(i)*[ones(1,KK) zeros(1,M-KK)];
end

figure(1);
index = length(a)+1;
t1 = Tb:Tb:(index-1)*Tb;
t2 = 0:Tb:(index-1)*Tb;
t3 = 0:ts:index*Tb-ts;
subplot(6, 1, 1);
stem(t1,a);grid
axis([0 index*Tb -1.2 1.2]);
title('输入信息码');ylabel('幅度');xlabel('t');

subplot(6, 1, 2);
stem(t2, b);grid
axis([0 index*Tb -1.2 1.2]);
title('传号差分码');ylabel('幅度');xlabel('t');

subplot(6, 1, 3);
plot(t3, b_s_NRZ(1:index*M));grid
axis([0 index*Tb -1.2 1.2]);
title('单极性不归零波形');ylabel('幅度');xlabel('t');

subplot(6, 1, 4);
plot(t3, b_s_RZ(1:index*M));grid
axis([0 index*Tb -1.2 1.2]);
title('单极性归零波形');ylabel('幅度');xlabel('t');

subplot(6, 1, 5);
plot(t3, b_d_NRZ(1:index*M));grid
axis([0 index*Tb -1.2 1.2]);
title('双极性不归零波形');ylabel('幅度');xlabel('t');

subplot(6, 1, 6);
plot(t3, b_d_RZ(1:index*M));grid
axis([0 index*Tb -1.2 1.2]);
title('双极性归零波形');ylabel('幅度');xlabel('t');

figure(2)
subplot(4, 1, 1);
stem(t1, aa);grid
axis([0 index*Tb -1.2 1.2]);
title('输入信息码');ylabel('幅度');xlabel('t');

subplot(4, 1, 2);
stem(t2, bb);grid
axis([0 index*Tb -1.2 1.2]);
title('传号差分码');ylabel('幅度');xlabel('t');

subplot(4, 1, 3);
plot(t3, bb_NRZ(1:index*M));grid
axis([0 index*Tb -1.2 1.2]);
title('双极性不归零波形');ylabel('幅度');xlabel('t');

subplot(4, 1, 4);
plot(t3, bb_RZ(1:index*M));grid
axis([0 index*Tb -1.2 1.2]);
title('双极性归零波形');ylabel('幅度');xlabel('t');
实验效果

单极性矩形波数字基带信号的功率谱仿真

实验代码
clc,clear;
% 参数设置
Tb = 1;
M = 20;
ts = Tb/M;
K = 0.5;
KK = floor(M*K);
N = 1000;
fs = 1/ts;
df = fs/N/M;
f = 0:df:(fs-df);

% 数值计算
P_lilun_NRZ = (sinc(f*Tb)).^2*Tb/4;
P_lilun_RZ = (sinc(f*Tb/2)).^2*Tb/16;

% 数据仿真
s = randi(1, N);

% 形成指定波形的基带信号
for i = 1:length(s)
    s_NRZ((i-1)*M+1:i*M) = s(i)*ones(1,M);
    s_RZ((i-1)*M+1:i*M) = s(i)*[ones(1,KK) zeros(1,M-KK)];
end
P_NRZ = (abs(fft(s_NRZ))/fs).^2;
P_RZ = (abs(fft(s_RZ))/fs).^2;

% 画图
index = length(s);
subplot(4, 1, 1);
semilogy(f, P_lilun_NRZ./max(P_lilun_NRZ));
axis([-0.1 6/Tb 10^-15 1.1]);
ylabel('归一化对数幅度');xlabel('f/Hz');title('不归零波形连续谱数值结果');grid

subplot(4, 1, 2);
semilogy(f, P_NRZ./max(P_NRZ));
axis([-0.1 6/Tb 10^-15 1.1]);
ylabel('归一化对数幅度');xlabel('f/Hz');title('不归零波形连续谱仿真结果');grid

subplot(4, 1, 3);
semilogy(f, P_lilun_RZ./max(P_lilun_RZ));
axis([-0.1 6/Tb 10^-15 1.1]);
ylabel('归一化对数幅度');xlabel('f/Hz');title('归零波形连续谱数值结果');grid


subplot(4, 1, 4);
semilogy(f, P_RZ./max(P_RZ));
axis([-0.1 6/Tb 10^-15 1.1]);
ylabel('归一化对数幅度');xlabel('f/Hz');title('归零波形连续谱仿真结果');grid
实验效果

双极性矩形波数字基带信号的功率谱仿真

实验代码
clc,clear;
% 参数设置
Tb = 1;
M = 20;
ts = Tb/M;
K = 0.5;
KK = floor(M*K);
N = 1000;
fs = 1/ts;
df = fs/N/M;
f = 0:df:(fs-df);

% 数值计算
P_lilun_NRZ = (sinc(f*Tb)).^2*Tb;
P_lilun_RZ = (sinc(f*Tb/2)).^2*Tb/4;

% 数据仿真
s = randi(1, N);

% 形成指定波形的基带信号
for i = 1:length(s)
    s_NRZ((i-1)*M+1:i*M) = (-1).^(s(i)-1)*ones(1,M);
    s_RZ((i-1)*M+1:i*M) = (-1).^(s(i)-1)*[ones(1,KK) zeros(1,M-KK)];
end
P_NRZ = (abs(fft(s_NRZ))/fs).^2;
P_RZ = (abs(fft(s_RZ))/fs).^2;

% 画图
index = length(s);
subplot(4, 1, 1);
semilogy(f, P_lilun_NRZ./max(P_lilun_NRZ));
axis([-0.1 6/Tb 10^-15 1.1]);
ylabel('归一化对数幅度');xlabel('f/Hz');title('不归零波形连续谱数值结果');grid

subplot(4, 1, 2);
semilogy(f, P_NRZ./max(P_NRZ));
axis([-0.1 6/Tb 10^-15 1.1]);
ylabel('归一化对数幅度');xlabel('f/Hz');title('不归零波形连续谱仿真结果');grid

subplot(4, 1, 3);
semilogy(f, P_lilun_RZ./max(P_lilun_RZ));
axis([-0.1 6/Tb 10^-15 1.1]);
ylabel('归一化对数幅度');xlabel('f/Hz');title('归零波形连续谱数值结果');grid

subplot(4, 1, 4);
semilogy(f, P_RZ./max(P_RZ));
axis([-0.1 6/Tb 10^-15 1.1]);
ylabel('归一化对数幅度');xlabel('f/Hz');title('归零波形连续谱仿真结果');grid
实验效果

计算零均值,单位方差的高斯随机分布的量化电平

实验代码
clc,clear;
m = 0;
sigma = 1;
a = [-10, -5, -4, -2, 0, 1, 3, 5, 10];

for i = 1:length(a)-1
    dx = (a(i+1) - a(i))/1000;
    x = a(i):dx:a(i+1);
    fx = normal_pdf(x, m, sigma);
    y1 = sum(fx)*dx;
    y2 = sum(x.*fx)*dx;
    y(i) = y2/y1;
end
y
实验效果

计算零均值,单位方差的高斯随机分布的量化电平量化器的量化噪声功率

实验代码
clc,clear;
m = 0;
sigma = 1;
a = [-10, -5, -4, -2, 0, 1, 3, 5, 10];

for i = 1:length(a)-1
    dx = (a(i+1) - a(i))/1000;
    x = a(i):dx:a(i+1);
    fx = normal_pdf(x, m, sigma);
    y1 = sum(fx)*dx;
    y2 = sum(x.*fx)*dx;
    y(i) = y2/y1;
    k = x - y(i);
    N(i) = sum(k.^2.*fx)*dx;
end
Nq = sum(N)
实验效果

计算均匀量化器对高斯分布信号的量化噪声功率

实验代码
clc,clear;
m = 0;
sigma = 4;
a = [-10, -5, -4, -2, 0, 1, 3, 5, 10];

for i = 1:length(a)-1
    dx = (a(i+1) - a(i))/1000;
    x = a(i)+dx:dx:a(i+1);
    fx = normal_pdf(x, m, sigma);
    y1 = sum(fx)*dx;
    y2 = sum(x.*fx)*dx;
    y(i) = y2/y1;
    k = x - y(i);
    N(i) = sum(k.^2.*fx)*dx;
end
y
Nq = sum(N)
实验效果

计算量化电平位于区域中点的均匀量化器对高斯分布信号的量化噪声功率

实验代码
clc,clear;
m = 0;
sigma = 4;
a = [-10, -5, -4, -2, 0, 1, 3, 5, 10];
y = [-7.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 7.5];

for i = 1:length(a)-1
    dx = (a(i+1) - a(i))/1000;
    x = a(i)+dx:dx:a(i+1);
    fx = normal_pdf(x, m, sigma);
    k = x - y(i);
    N(i) = sum(k.^2.*fx)*dx;
end
Nq = sum(N)
实验效果

高斯随机变化序列的均匀PCM量化

实验代码
clc,clear;
% 产生高斯随机变化序列
x = randn(1, 500);
% 量化与编码
n = 64;
[xq, code_n, delta] = unipcm(x, n);
% 计算量化信噪比
Nq = delta^2/12;
Sq = (norm(xq))^2/length(xq);
sqnr = 10*log10(Sq/Nq);
index = 5;
x(1:index)
xq(1:index)
code_n(1:index*log2(n))
% 画图
subplot(2, 1, 1)
t = 1:1:length(x);
plot(t, x-xq, 'r');hold on
xlabel('样值序列');ylabel('幅度');title('量化误差');
subplot(2, 1, 2)
[y, I] = sort(x);
plot(y, xq(I));grid
xlabel('输入序列值');ylabel('幅度');title('以输入序列为函数的量化输出信号');
实验效果

A律13折线仿真

实验代码
clc,clear;
% 对x轴部分
x(1) = 0;
x(9) = 1;
for i = 7:-1:1
    x(9-i) = 1/(2^i);
end

% 对y轴部分
y = 0:1/8:1;

% A律特性
A = 87.6;
dt = 1/1000;
xx = 0:dt:1;
for i = 1:length(xx)
    if xx(i) <= 1/A;
        yy(i) = A*xx(i)/(1+log(A));
    else
        yy(i) = (1+log(A*xx(i)))/(1+log(A));
    end
end

% 画图
plot(xx, yy);hold on
plot(x, y, 'ro--');hold on
xlabel('归一化输入');ylabel('归一化输出');title('A律压缩特性');grid
legend('A律压缩特性','A律13折线')
实验效果

两路径无线信道模型仿真

实验代码
clc,clear;
% 参数设置
p = [0.99 0.9];
N = 1000;
d = 5;             %延时样本数
for i = 1:length(p)
    A = [1 -2*p(i) p(i)^2];        %滤波器参数
    B = (1-p(i))^2;
    white_noise_seq1 = randn(1,N);
    white_noise_seq2 = randn(1,N);
    R1(i, :) = filter(B, A, white_noise_seq1);
    R2(i, :) = filter(B, A, white_noise_seq2);
    h(i, :) = R1(i,d+1:N)+R2(i,1:N-d);
end
subplot(2, 1, 1);
plot(1:N,R1(1, :),'-',1:N,R2(1, :),'--',1:N-d,h(1, :),':');
legend('径1', '径2', '合成信道')
xlabel('n');ylabel('幅度');title('p=0.99')
subplot(2, 1, 2);
plot(1:N,R1(2, :),'-',1:N,R2(2, :),'--',1:N-d,h(2, :),':');
legend('径1', '径2', '合成信道')
xlabel('n');ylabel('幅度');title('p=0.9')
实验效果

产生受遮蔽的Rician衰落分布的信道系数

实验代码
clc,clear;
% 产生对数分布均值的复高斯随机变量
N = 10^5;     % 样值点数
sigma_a = 0.252;    % 散射信号分量方差
x = sqrt(sigma_a/2)*randn(2, N);
sigma_z = 0.161^2;   % 直射信号分量对数的方差
u_z = -0.115;        % 直射信号分量对数的均值
zz = sqrt(sigma_z)*randn(1, N);
z = exp(zz +u_z);
xx = z + x(1, :) + sqrt(-1)*x(2, :);
k = 0:0.01:5;
pdf_sim = c_pdf(k, abs(xx));
% 理论值
kk = 0:0.01:5;
dz = 0.001;
zk = dz:dz:10;
for i = 1:length(kk)
    xa = 2*kk(i)/(sqrt(2*pi*sigma_z)*sigma_a);
    xb = exp(-(kk(i)^2+zk.^2)/sigma_a-(log(zk)-u_z).^2/(2*sigma_z));
    xc = besseli(0, 2*kk(i)*zk/sigma_a)*dz./zk;
    xd = sum(xb.*xc);
    pdf_lilun(i) = xa*xd;
end
% 画图
plot(kk, pdf_lilun, 'r');hold on
plot(k, pdf_sim, '-');grid
xlabel('R');ylabel('PDF');legend('理论值', '仿真值');
title('中度遮蔽的受遮蔽Rician衰落分布')
实验效果

画出低通白噪声的自相关函数曲线

实验代码
clc,clear;
bw = 1;
fs = 100;
ts = 1/fs;
tao = -5:ts:5;
N0 = 1;
B = 1;
Rx = N0*B*sinc(2*B*tao);
figure(1);
plot(tao, Rx);grid on;
xlabel('tao/s');ylabel('自相关');title('低通白噪声的自相关函数');

% 求自相关函数和功率谱
N = 10000;
M = 500;
Rxav = zeros(1,M+1);
Ryav = zeros(1,M+1);
Sxav = zeros(1,M+1);
Syav = zeros(1,M+1);
% 设计带宽为BW = 1Hz的低通滤波器
bw = 1;
a = 1;
b1 = bw/(fs/2);
b = fir1(1000, b1);
% 求自相关函数及功率谱密度
for i =1:10
    detla2 = N0*(fs/2);
    se_orignal = sqrt(detla2)*randn(1, N);    %产生1000个单边功率谱密度为N0=1W/Hz的白噪声序列
    low_out = filter(b, a, se_orignal);
    % 求se_orignal序列的自相关函数
    Rx = Rx_est(se_orignal, M);
    Ry = Rx_est(low_out, M);
    Sx = fftshift(abs(fft(Rx)/fs));
    Sy = fftshift(abs(fft(Ry)/fs));
    Rxav = Rxav + Rx;
    Ryav = Ryav + Ry;
    Sxav = Sxav + Sx;
    Syav = Syav + Sy;
end
Rxav = Rxav/10;
weight = Rxav(1);
Rxav_1 = Rxav/weight;    % 归一化
Ryav = Ryav/10;
weight = Ryav(1);
Ryav_1 = Ryav/weight;    % 归一化
Sxav = Sxav/10;
Syav = Syav/10;
figure(2);
t = (-M: M)*ts;
plot(t, [Rxav_1(M:-1:1), Rxav_1]);grid;  %画出归一化白噪声自相关函数
xlabel('tao/s');ylabel('归一化自相关');title('白噪声的仿真自相关函数');
figure(3);
plot(t, [Ryav_1(M:-1:1), Ryav_1]);grid;  %画出归一化低通白噪声自相关函数
xlabel('tao/s');ylabel('归一化自相关');title('低通白噪声的仿真自相关函数');
figure(4);
df = fs/length(Sxav);
f = (0:length(Sxav)-1)*df-fs/2;
plot(f, Sxav);     % 画出白噪声的频谱(其平方为功率谱)
xlabel('f/Hz');ylabel('幅度');title('白噪声仿真频谱')
figure(5);
plot(f, Syav);     % 画出低通白噪声的频谱(其平方为功率谱)
xlabel('f/Hz');ylabel('幅度');title('低通白噪声仿真频谱')
实验效果

解说题3.5 DSB-AM解调

实验代码
clc,clear;
echo on
t0=.15;                                 % 信号持续时间
ts=1/1500;                              % 采样间隔
fc=250;                                 % 载波频率
fs=1/ts;                                % 采样频率
t=[0:ts:t0];                            % 时间矢量
df=0.3;                                 % 所需频率分辨率
% message signal
m=[ones(1,t0/(3*ts)),-2*ones(1,t0/(3*ts)),zeros(1,t0/(3*ts)+1)];
c=cos(2*pi*fc.*t);                      % 载波信号
u=m.*c;                                 % 调制信号
y=u.*c;                                 % 混合信号
[M,m,df1]=fftseq(m,ts,df);              % FFT
M=M/fs;                                 % 缩放比例
[U,u,df1]=fftseq(u,ts,df);              % FFT
U=U/fs;                                 % 缩放比例
[Y,y,df1]=fftseq(y,ts,df);              % FFT
Y=Y/fs;                                 % 缩放比例
f_cutoff=150;                           % 滤波器的截止频率
n_cutoff=floor(150/df1);                % 设计过滤器。
f=[0:df1:df1*(length(y)-1)]-fs/2;
H=zeros(size(f));
H(1:n_cutoff)=2*ones(1,n_cutoff);
H(length(f)-n_cutoff+1:length(f))=2*ones(1,n_cutoff);
DEM=H.*Y;                               % 滤波器输出的频谱
dem=real(ifft(DEM))*fs;                 % 滤波器输出
subplot(4,3,1)
plot(f,fftshift(abs(M)))
title('消息信号频谱')
xlabel('频率')
subplot(4,3,2)
plot(f,fftshift(abs(U)))
title('已调信号的频谱')
xlabel('频率')
subplot(4,3,3)
plot(f,fftshift(abs(Y)))
title('混合输出频谱')
xlabel('频率')
subplot(4,3,4)
plot(f,fftshift(abs(Y)))
title('混合输出频谱')
xlabel('频率')
subplot(4,3,5)
plot(f,fftshift(abs(H)))
title('低通滤波器频谱')
xlabel('频率')
subplot(4,3,6)
plot(f,fftshift(abs(DEM)))
title('滤波后的频谱')
xlabel('频率')
subplot(4,3,7)
plot(f,fftshift(abs(M)))
title('消息信号的频谱')
xlabel('频率')
subplot(4,3,9)
plot(f,fftshift(abs(DEM)))
title('滤波后的频谱')
xlabel('频率')
subplot(4,3,10)
plot(t,m(1:length(t)))
title('消息信号')
xlabel('时间')
subplot(4,3,12)
plot(t,dem(1:length(t)))
title('恢复后的信号')
xlabel('时间')
实验效果

解说题3.11 频率调制

实验代码
clc,clear;
echo on
clear;
t0=.2;                               	% signal duration
ts=0.001;                            	% sampling interval
fc=250;                              	% carrier frequency
snr=20;                              	% SNR in dB (logarithmic)
fs=1/ts;                             	% sampling frequency
df=0.3;                              	% required freq. resolution
t=[-t0/2:ts:t0/2];                   	% time vector
kf=100;                              	% deviation constant
df=0.25;                             	% required frequency resolution
m=sinc(100*t);                       	% 消息信号
int_m(1)=0;
for i=1:length(t)-1                  	% integral of m
    int_m(i+1)=int_m(i)+m(i)*ts;
    echo off ;
end
echo on ;
[M,m,df1]=fftseq(m,ts,df);           	% 傅里叶变换
M=M/fs;                              	% scaling
f=[0:df1:df1*(length(m)-1)]-fs/2;    	% frequency vector
u=cos(2*pi*fc*t+2*pi*kf*int_m);      	% 已调信号
[U,u,df1]=fftseq(u,ts,df);           	% Fourier transform
U=U/fs;                              	% scaling
[v,phase]=env_phas(u,ts,250);        	% 解调,找到u的相位
phi=unwrap(phase);                   	% 还原原始的相位
dem=(1/(2*pi*kf))*(diff(phi)/ts);    	% 解调输出, differentiate and scale phase
% pause  % Press any key to see a plot of the message and the modulated signal.
subplot(3,2,1)
plot(t,m(1:length(t)))
xlabel('时间')
title('消息信号')
subplot(3,2,2)
plot(t,u(1:length(t)))
xlabel('时间')
title('已调信号')
% pause   % Press any key to see plots of the magnitude of the message and the
% modulated signal in the frequency domain.
subplot(3,2,3)
plot(f,abs(fftshift(M)))
xlabel('频率')
title('消息信号的频谱')
subplot(3,2,4)
plot(f,abs(fftshift(U)))
title('已调信号的频谱')
xlabel('频率')
% pause	% Press any key to see plots of the message and the demodulator output with no
% noise.
subplot(3,2,5)
plot(t,m(1:length(t)))
xlabel('时间')
title('消息信号')
subplot(3,2,6)
plot(t,dem(1:length(t)))
xlabel('时间')
title('解调信号')
实验效果

解说题7.9 QAM信号的解调

实验代码
clc,clear;
M = 8;
Es = 1;                                 
T = 1;
Ts = 100/T;
fc = 30/T;
t = 0:T/100:T;
l_t = length(t);
A_mc = 1/sqrt(Es);                      % 信号幅度
A_ms = -1/sqrt(Es);                     % 信号幅度
g_T = sqrt(2/T)*ones(1,l_t);
phi = 2*pi*rand;
si_1 = g_T.*cos(2*pi*fc*t + phi);
si_2 = g_T.*sin(2*pi*fc*t + phi);
var = [ 0 0.05 0.5];                    % 噪声方差矢量
for k = 1 : length(var)
    % 噪声成分的产生:
    n_c = sqrt(var(k))*randn(1,l_t);
    n_s = sqrt(var(k))*randn(1,l_t);
    noise = n_c.*cos(2*pi*fc+t) - n_s.*sin(2*pi*fc+t);
    % 接收到的信号
    r = A_mc*g_T.*cos(2*pi*fc*t+phi) + A_ms*g_T.*sin(2*pi*fc*t+phi) + noise;
    % 相关器输出:
    y_c = zeros(1,l_t);
    y_s = zeros(1,l_t);
    for i = 1:l_t
        y_c(i) = sum(r(1:i).*si_1(1:i));
        y_s(i) = sum(r(1:i).*si_2(1:i));
    end
    % 绘制结果:
    subplot(3,1,k)
    plot([0 1:length(y_c)-1],y_c,'.-')
    hold
    plot([0 1:length(y_s)-1],y_s)
    title(['\sigma^2 = ',num2str(var(k))])
    xlabel('n')
    axis auto
end
实验效果

解说题7.10 QAM信号的仿真

实验代码
clc,clear;echo on
SNRindB1=0:2:15;
SNRindB2=0:0.1:15;
M=16;
k=log2(M);
for i=1:length(SNRindB1)
    smld_err_prb(i)=cm_sm41(SNRindB1(i));	% 模拟错误率
    echo off;
end
echo on ;
for i=1:length(SNRindB2)
    SNR=exp(SNRindB2(i)*log(10)/10);    	% 信噪比
    % 理论符号错误率
    theo_err_prb(i)=4*Qfunct(sqrt(3*k*SNR/(M-1)));
    echo off ;
end
echo on ;
% 绘图
semilogy(SNRindB1,smld_err_prb,'*');
hold
semilogy(SNRindB2,theo_err_prb);
实验效果

解说题8.2 OFDM信号的产生

实验代码
clc,clear;echo on
K=10;N=2*K;T=100;
a=rand(1,36);
a=sign(a-0.5);
b=reshape(a,9,4);
% 生成16QAM点
XXX=2*b(:,1)+b(:,2)+1i*(2*b(:,3)+b(:,4));
XX=XXX';
X=[0 XX 0 conj(XX(9:-1:1))];
xt=zeros(1,101);
for t=0:100
    for k=0:N-1
        xt(1,t+1)=xt(1,t+1)+1/sqrt(N)*X(k+1)*exp(1i*2*pi*k*t/T);
        echo off
    end
end
echo on
xn=zeros(1,N);
for n=0:N-1
    for k=0:N-1
        xn(n+1)=xn(n+1)+1/sqrt(N)*X(k+1)*exp(1i*2*pi*n*k/N);
        echo off
    end
end
echo on
plot(0:100,abs(xt))
% 检查xn与x(t)样本之间的差异
for n=0:N-1
    d(n+1)=xt(T/N*n+1)-xn(1+n);
    echo off
end
echo on
e=norm(d);
Y=zeros(1,10);
for k=1:9
    for n=0:N-1
        Y(1,k+1)=Y(1,k+1)+1/sqrt(N)*xn(n+1)*exp(-1i*2*pi*k*n/N);
        echo off
    end
end
echo on
dd=Y(1:10)-X(1:10);
ee=norm(dd);
实验效果

解说题8.5 OFDM信号的频谱

实验代码
clc,clear;
T = 1;
k = 0 : 5;
f_k = k/T;
f = 0 : 0.01*4/T : 4/T;
U_k_abs = zeros(length(k),length(f));
for i = 1 : length(k)
    U_k_abs(i,:) = abs(sqrt(T/2)*(sinc((f-f_k(i))*T) + sinc((f+f_k(i))*T)));
    U_k_norm(i,:) = U_k_abs(i,:)/max(U_k_abs(i,:));
end
U_k_dB = 10*log10(U_k_norm);
plot(f,U_k_dB(1,:),'black',f,U_k_dB(2,:),'black',f,U_k_dB(3,:),'black',...
    f,U_k_dB(4,:),'black',f,U_k_dB(5,:),'black',f,U_k_dB(6,:),'black')
axis([min(f) max(f) -180 20])
xlabel('f')
ylabel('|U_k(f)| (dB)')
实验效果

全部子函数

实验代码
%% 高斯分布函数

function[y] = normal_pdf(x, m, sigma)
% 产生均值为m,方差为sigma的高斯分布函数
% x取值区间为[a b]

y = exp(-(x-m).^2/(2*sigma))/sqrt(2*pi*sigma);



%% 均与量化与自然编码子程序,x -- 模拟信源, n -- 量化级数
function[x_q, code, delat] = unipcm(x, n)
M = log2(n);

x_width = max(x) - min(x);                                 % 量化范围的大小 
delat = x_width/n;                                         % 量化间隔
xx =  min(x):delat:max(x);                                 % 量化分层电平值   
q = min(x) + delat/2:delat:max(x);                         % 量化电平值

for i = 1:length(x)
    if x(i) <= xx(2)                                       % 小于等于第2个量化分层的均属于第1量化级
        x_q(i) = q(1);
        index(i) = 1;
    end
    k = 2;
    while k < n+1
        if(x(i) > xx(k)) & (x(i)<=xx(k+1))                 % xx(k) < x(i) <= xx(k+1),位于分层电平
            index(i) = k;
            x_q(i) = q(k);
            k = n+1;
        end
        k = k+1;
    end
    code(i*M:-1:(i-1)*M+1) = de2bi(index(i)-1, M);         % 二进制编码
end



function pdf = c_pdf(k, data)
    [counts, binCenters] = hist(data, k);
    pdf = counts / (sum(counts) * (binCenters(2) - binCenters(1)));
end



function Rx = Rx_est(x, M)
    % 计算自相关函数
    N = length(x);
    Rx = zeros(1, M+1);
    for m = 0:M
        Rx(m+1) = sum(x(1:N-m) .* x(m+1:N)) / N;
    end
end



function xl=loweq(x,ts,f0)
%       xl=loweq(x,ts,f0)
%LOWEQ      returns the lowpass equivalent of the signal x
%       f0 is the center frequency. 
%       ts is the sampling interval.
%   
t=[0:ts:ts*(length(x)-1)];  
z=hilbert(x);
xl=z.*exp(-j*2*pi*f0*t);



function [M,m,df]=fftseq(m,ts,df) 
%       [M,m,df]=fftseq(m,ts,df)
%       [M,m,df]=fftseq(m,ts)
%FFTSEQ     generates M, the FFT of the sequence m.
%       The sequence is zero padded to meet the required frequency resolution df.
%       ts is the sampling interval. The output df is the final frequency resolution.
%       Output m is the zero padded version of input m. M is the FFT.
fs=1/ts;
if nargin == 2
  n1=0;
else
  n1=fs/df;
end
n2=length(m);
n=2^(max(nextpow2(n1),nextpow2(n2)));
M=fft(m,n);
m=[m,zeros(1,n-n2)];
df=fs/n;



function [v,phi]=env_phas(x,ts,f0)
%		[v,phi]=env_phas(x,ts,f0)
%		      v=env_phas(x,ts,f0)
%ENV_PHAS   	Returns the envelope and the phase of the bandpass signal x.
%		f0 is the center frequency. 
%		ts is the sampling interval.
%
if nargout == 2
  z=loweq(x,ts,f0);
  phi=angle(z);
end
v=abs(hilbert(x));