function n_t = bandlimit_noise(N,fs,fc,B,n_pow)
% 产生特定功率的带限高斯白噪声
%  输入： N  - 产生的噪声的长度
%        fs - 仿真采样率
%        fc - 带限噪声中心频率
%        B  - 带限噪声的带宽
%        n_pow - 带限噪声的功率
%  输出： n_t - 带限噪声%带通阶数
filter_b_L = 4000;
%带通设计
filter_b = fir1(4000,[(fc - B/2) / (fs/2) , (fc + B/2) / (fs/2)]);
%全带宽功率
n_fullband_pow = n_pow/(B/(fs/2));
%产生噪声信号
noise_t = sqrt(n_fullband_pow) * randn( N + filter_b_L,1);
%信号经过带通滤波器
n_t = filter(filter_b,1,noise_t);
%滤除前面的零状态响应
n_t = n_t(end-N+1:end);end

% 产生特定功率的带限高斯白噪声
%  输入： N  - 产生的噪声的长度
%        fs - 仿真采样率
%        fc - 带限噪声中心频率
%        B  - 带限噪声的带宽
%        n_pow - 带限噪声的功率
%  输出： n_t - 带限噪声clc,clear;N = 5e5;
fs = 20e3;
fc = 3e3;
B = 400;
n_pow = 0.1;t = (1:N)/fs;
noise = bandlimit_noise(N,fs,fc,B,n_pow);plot(t,noise);var(noise)

function snr_est = sine_snr_est(y)
% 正弦波信噪比估计
%        输入信号序列；输出为信号信噪比(单位dB)N = floor(length(y)/2) * 2;y = y(1:N);y = y.*blackman(N);f_n = abs(fft(y));f_n = f_n(1:N/2);f_pow = abs(f_n).^2;[s_pow_est idx]= max(f_pow);s_pow_idx = [idx-4:idx+4];s_pow_est = sum(f_pow(s_pow_idx));n_pow_est = sum(f_pow) - s_pow_est;snr_est = db(s_pow_est/n_pow_est,'power');end

% 根据参数产生AM信号
%  输入： N  - 产生的信号长度
%        fs - 仿真采样率
%        fm - 调制信号频率
%        Am - 调制信号幅度
%        yita - 调幅系数
%        fc - 载波频率
%        Ac - 载波幅度
%  输出：s_t - 已调信号
%        s_pow - 已调信号功率clc,clear;%噪声（第一题的）
N = 5e4;
fs = 20e3;
fc = 3e3;
B = 400;
n_pow = 0.1;%正弦信号
t = ((0:N-1) / fs)';
m_t = sin( 2*pi*fc*t );
m_t_pow = sum(m_t.^2) / N;%产生高斯带通噪声
SNR_in = -10:40;
n_pow = m_t_pow./(10.^(SNR_in/10));
for i=1:length(n_pow)n_t = bandlimit_noise(N,fs,fc,B,n_pow(i));sign_in = m_t + n_t;SNR_out(i) = sine_snr_est(sign_in);
endplot(SNR_in,SNR_out);
xlabel('实际SNR(Si/Ni)');
ylabel('检测SNR(So/No)');

function snr_db_o = am_receiver(snr_db_i,flag)
% 根据输入信噪比snr_db_i，通过仿真确定AM解调器输出信噪比% 仿真参数fs = 20e3;N = 1000000;% 载波参数fc = 2e3;Ac = 1;% 调制信号参数Am = 1;bw = 300;fm = bw/2;yita = 1;% 产生调制信号[s_t,s_pow] = am_modulate(N,fs,fm,Am,yita,fc,Ac);% 加入噪声n_pow = s_pow/(10^(snr_db_i/10));n_t = bandlimit_noise(N,fs,fc,bw,n_pow);s_t = s_t + n_t;% 解调if flag==1r_demod_t = am_demodulate_xianggan(fs,fc,bw,s_t);elser_demod_t = am_demodulate_feixianggan(fs,fc,bw,s_t);end% 输出信号SNR估计snr_db_o = sine_snr_est(r_demod_t);end

function [s_t,s_pow] = am_modulate(N,fs,fm,Am,yita,fc,Ac)
% 根据参数产生AM信号
%  输入： N  - 产生的信号长度
%        fs - 仿真采样率
%        fm - 调制信号频率
%        Am - 调制信号幅度
%        yita - 调幅系数
%        fc - 载波频率
%        Ac - 载波幅度
%  输出：s_t - 已调信号
%        s_pow - 已调信号功率% 时间t,作为横坐标
t = (0:N-1)'/fs;
% 产生调制信号
m_t = Am * sin(2*pi*fm*t);
% 产生载波信号
c_t = Ac * cos(2*pi*fc*t);
% 调制
A0 = Am / yita;
s_t = (m_t + A0) .* c_t; s_pow = sum(s_t.^2)/N;end

function r_demod_t = am_demodulate_xianggan(fs,fc,bw,s_t)
% AM信号解调，相干
%  输入： fs - 仿真采样率
%        fc - 载波频率
%        bw - 基带带宽
%        s_t - 已调信号
%  输出： r_demod_t - 解调后的信号r_t = s_t;
N = length(s_t);
t = (0:N-1)'/fs;
c_t = cos(2*pi*fc*t); % 本地振荡器波形
mix_t = r_t .* c_t; % 混频% 滤波
b = fir1(20, bw/(fs/2)); % 设计截止频率为fcFIR低通滤波器
filtered_t = filter(b,1,mix_t); % 去除DC分量
r_demod_t = filtered_t - sum(filtered_t)/length(filtered_t);end

function r_demod_t = am_demodulate_feixianggan(fs,fc,bw,s_t)
% AM信号解调，非相干
%  输入： fs - 仿真采样率
%        fc - 载波频率
%        bw - 基带带宽
%        s_t - 已调信号
%  输出： r_demod_t - 解调后的信号%半波整流
s_t = s_t .* (s_t > 0);% 滤波
b = fir1(80, bw/(fs/2));           % 设计截止频率为fc/2的FIR低通滤波器
filtered_t = filter(b,1,s_t);       %去直流
r_demod_t = filtered_t - sum(filtered_t)/length(filtered_t);end

% 对AM接收机进行性能评估
clc,clear;% 输入SNR从-5dB变化到25dBsnr_db_i = -5:1:25;snr_db_o = zeros(1,length(snr_db_i));for k = 1:length(snr_db_i)snr = snr_db_i(k);snr_db_o(k) = am_receiver(snr,1);end% 绘制图形figure; plot(snr_db_i,snr_db_o); title('AM 相干解调器的“解调输出SNR(So/No) - 输入SNR(Si/Ni)”曲线'); xlabel('输入SNR(Si/Ni)'); ylabel('输出SNR(So/No)');% 输入SNR从-5dB变化到25dBsnr_db_i = -5:1:25;snr_db_o = zeros(1,length(snr_db_i));for k = 1:length(snr_db_i)snr = snr_db_i(k);snr_db_o(k) = am_receiver(snr,2);end% 绘制图形figure; plot(snr_db_i,snr_db_o); title('AM非相干解调器的“解调输出SNR(So/No) - 输入SNR(Si/Ni)”曲线'); xlabel('输入SNR(Si/Ni)'); ylabel('输出SNR(So/No)');

function snr_db_o = fm_receiver(snr_db_i,mf)
% 根据输入信噪比snr_db_i，通过仿真确定AM解调器输出信噪比% 仿真参数fs = 20e3;N = 1000000;% 载波参数fc = 2e3;Ac = 1;% 调制信号参数Am = 1;bw = 300;fm = bw/2;% 产生调制信号[s_t,s_pow] = fm_modulate(N,fs,fm,Am,mf,fc,Ac);% 加入噪声n_pow = s_pow/(10^(snr_db_i/10));n_t = bandlimit_noise(N,fs,fc,bw,n_pow);s_t = s_t + n_t;% 解调r_demod_t = demod(s_t,fc,fs,'fm');% 输出信号SNR估计snr_db_o = sine_snr_est(r_demod_t);end

function [s_t,s_pow] = fm_modulate(N,fs,fm,Am,mf,fc,Ac)
% 根据参数产生FM信号
% fm      调制信号频率(Hz)
% Am      调制信号幅度
% fc      载波频率(Hz)
% Ac      载波幅度
% mf      调频指数
% N       数据点数
% fs      采样频率kf = mf * 2 * pi * fm / Am;t = (0:N-1)'/fs; % 时间t
m_t = Am*sin(2*pi*fm*t); % 调制信号
phi_t = kf*cumsum(m_t)/fs; % 相位积分
s_t = cos(2*pi*fc*t + phi_t); % 已调信号s_pow = sum(s_t.^2) / N;%功率计算end

% 对FM接收机进行性能评估
clc,clear;% 输入SNR从-5dB变化到25dBsnr_db_i = 5:1:20;snr_db_o = zeros(1,length(snr_db_i));figure;hold on;for k = 1:length(snr_db_i)snr = snr_db_i(k);snr_db_o(k) = fm_receiver(snr,0.5);end% 绘制图形plot(snr_db_i,snr_db_o,'r'); for k = 1:length(snr_db_i)snr = snr_db_i(k);snr_db_o(k) = fm_receiver(snr,1);end% 绘制图形plot(snr_db_i,snr_db_o,'g'); for k = 1:length(snr_db_i)snr = snr_db_i(k);snr_db_o(k) = fm_receiver(snr,2);end% 绘制图形plot(snr_db_i,snr_db_o,'b'); for k = 1:length(snr_db_i)snr = snr_db_i(k);snr_db_o(k) = fm_receiver(snr,5);end% 绘制图形plot(snr_db_i,snr_db_o,'k'); for k = 1:length(snr_db_i)snr = snr_db_i(k);snr_db_o(k) = fm_receiver(snr,7);end% 绘制图形plot(snr_db_i,snr_db_o,'y');legend('mf=0.5','mf=1','mf=2','mf=5','mf=7');title('FM 相干解调器的“解调输出SNR(So/No) - 输入SNR(Si/Ni)”曲线'); xlabel('输入SNR(Si/Ni)'); ylabel('输出SNR(So/No)');