%以下代码功能：1.读取音频文件并对音频文件进行低通滤波，截止频率9000Hz；
% 2.生成滤波后的音频文件试听；
% 3.对滤波前后音频文件时频域进行分析；
% 4.对滤波后音频进行预加重并在时频域进行分析
% 5.将两个音频文件信号重采样为3MHz
% 6.将音频信号进行fm调制
% 7.调制信号AWGN信道
% 8.调用fmdemod函数对调制信号进行解调
% 9.对解调后的信号进行重采样为44.1KHz
% 10.对恢复信号进行分析
% 11.绘图分析
clc;clear;close all;
%设置文件路径
file_path='C:\Users';%需要设置为自己的路径
%设置信道信噪比
snr_set = 30;
%读取并分析音频时频域
[x,Fs]=audioread([file_path,'\wheniwasyoung.mp3']);
x = x(:,1);
x = x'; %转置，这里把列向量变成行向量
x =x/max(abs(x)); %归一化
%准备变采样前的时间、频率变量，为绘图做准备
N = length(x); %求取抽样点数
t = (0:N-1)/Fs; %显示实际时间
f = Fs/N*(0:round(N/2)-1); %显示实际频点的一半，频域映射，转化为HZ
%----------------------------滤波------------------------------------------
% 设计一个低通FIR滤波器
n = 150; %滤波器的阶数
f_cutoff = 9000 / (Fs/2); %截止频率，由于fir1的规定，这里需要归一化
b = fir1(n, f_cutoff); %fir滤波器的分子系数
%分析滤波前信号的频谱
y= fft(x);
magnitude = abs(y(1:round(N/2)));
magnitude_db=20*log10(magnitude); %为画对数谱做准备
% 使用filter函数将滤波器应用到信号x上，并对滤波后的信号进行fft
x_filtered = filter(b, 1, x);
x_filtered = x_filtered/max(abs(x_filtered)); %滤波后再次归一化
%分析滤波后信号的频谱
y_filtered = fft(x_filtered);
magnitude_filtered = abs(y_filtered(1:round(N/2)));
magnitude_filtered_db=20*log10(magnitude_filtered); %为画对数谱做准备
%-----------------------预加重----------------------------------------------
%对信号进行预加重
x_preemph = filter([1 -0.9495],1,x_filtered);
x_preemph = x_preemph/max(abs(x_preemph)); %预加重后再次归一化
%分析预加重后信号频谱
y_preemph = fft(x_preemph);
magnitude_preemph = abs(y_preemph(1:round(N/2)));
magnitude_preemph_db=20*log10(magnitude_preemph); %为画对数谱做准备
%-----------------------重采样----------------------------------------------
%由于加入载波，对信号重采样
Fs_new = 3e6; % 重采样频率
x_resampled = resample(x_filtered, 68,1); %重采样原始音频信号 44.1KHz→3MHz
%注：解调后恢复采样率则是resample(解调信号,1,68) 3MHz→44.1KHz
N_resampled = length(x_resampled); %求重采样后信号点数
t_resampled = (0:N_resampled-1)/Fs_new; %重新设置的时间t
f_resampled = Fs_new/N_resampled*(0:round(N_resampled/2)-1);%显示实际频点的一半，频域映射，转化为HZ
%对预加重之后的信号重采样
x_preemph_resampled = resample(x_preemph, 68,1);
%-----------------------调制-----------------------------------------------
%fm调制参数设置
fc = 1e6; %设置载波频率
bandwidth = 9000; %设置基带带宽是9Khz
%根据预设参数计算kfm和βfm
beta = 9.659; %指定调频系数
kfm = bandwidth*beta; %计算频偏系数(音频信号已归一化，故最大幅度就是1)
%对信号进行直接调制
cal_t = cumtrapz(t_resampled,x_resampled); %对音频信号进行积分
phi =2*pi*fc*t_resampled +2*pi*kfm*cal_t;
s_fm = cos(phi); %生成调制信号
%对预加重之后的信号进行直接调制
cal_t_pre = cumtrapz(t_resampled,x_preemph_resampled); %对音频信号进行积分
phi_pre =2*pi*fc*t_resampled +2*pi*kfm*cal_t_pre;
s_fm_pre = cos(phi_pre); %生成调制信号
%--------------------------------AWGN信道----------------------------------
%对调制信号加噪声 并检验
s_fm_noise = awgn(s_fm,snr_set,'measured');
s_fm_pre_noise = awgn(s_fm_pre,snr_set,'measured');
%计算信噪比
SNR_in = snr_ana(s_fm,s_fm_noise);
SNR_in_pre = snr_ana(s_fm_pre,s_fm_pre_noise);
%-----------------------对调制信号进行分析及解调------------------------------------------
%---------------以下分析未加噪声、未使用预加重调制信号的频谱------------------
y_new = fft(s_fm); %fft
magnitude_sfm = abs(y_new(1:round(N_resampled/2))); %求幅度
f_new = Fs_new/N_resampled*(0:round(N_resampled/2)-1); %显示实际频点的一半，频域映射，转化为HZ
magnitude_sfm_db =20*log10(magnitude_sfm); %为画对数谱做准备
%解调
demodulated_signal = fmdemod(s_fm, fc, Fs_new, kfm); %调用自带解调函数，由于信号归一化，最大频偏就是kfm
demodulated_signal_resampled = resample(demodulated_signal, 1, 68);%重采样恢复原始音频信号
demodulated_signal_resampled = demodulated_signal_resampled / max(abs(demodulated_signal_resampled)); %归一化
%再计算解调后信号的频谱
y_de =fft(demodulated_signal_resampled); %对解调后的信号进行傅里叶变换
magnitude_de = abs(y_de(1:round(N/2))); %取幅度
magnitude_de_db=20*log10(magnitude_de); %为画对数谱做准备
%-------------------以下分析未加噪声、使用预加重调制信号的频谱----------------
y_new_pre = fft(s_fm_pre); %fft
magnitude_sfm_pre = abs(y_new_pre(1:round(N_resampled/2))); %求幅度
magnitude_sfm_pre_db =20*log10(magnitude_sfm_pre); %为画对数谱做准备
% 以下解调
demodulated_signal_pre_beforede = fmdemod(s_fm_pre, fc, Fs_new, kfm); %调用自带解调函数，由于信号归一化，最大频偏就是kfm
demodulated_signal_resampled_pre_beforede = resample(demodulated_signal_pre_beforede, 1, 68);%重采样恢复原始音频信号
demodulated_signal_resampled_pre_beforede = demodulated_signal_resampled_pre_beforede / max(abs(demodulated_signal_resampled_pre_beforede)); %归一化
%进行去加重
demodulated_signal_resampled_pre = filter(1, [1 -0.9495], demodulated_signal_resampled_pre_beforede);%去加重
demodulated_signal_resampled_pre = demodulated_signal_resampled_pre / max(abs(demodulated_signal_resampled_pre)); %归一化
% audiowrite([file_path,'\预加重但未去加重测试解调效果.wav'], demodulated_signal_resampled_pre_beforede, Fs); %解调的声音存储为文件
% audiowrite([file_path,'\预加重去加重测试解调效果.wav'], demodulated_signal_resampled_pre, Fs); %解调的声音存储为文件
%再计算解调后信号的频谱
%去加重前
y_de_pre_beforede =fft(demodulated_signal_resampled_pre_beforede);%对解调后的信号进行傅里叶变换
magnitude_de_pre_beforede = abs(y_de_pre_beforede(1:round(N/2))); %取幅度
magnitude_de_pre_beforede_db=20*log10(magnitude_de_pre_beforede);%为画对数谱做准备
%去加重后
y_de_pre =fft(demodulated_signal_resampled_pre);%对解调后的信号进行傅里叶变换
magnitude_de_pre = abs(y_de_pre(1:round(N/2))); %取幅度
magnitude_de_pre_db=20*log10(magnitude_de_pre);%为画对数谱做准备
%-------------------以下分析加噪声、未使用预加重调制信号的频谱----------------
y_new_noise = fft(s_fm_noise); %fft
magnitude_sfm_noise = abs(y_new_noise(1:round(N_resampled/2))); %求幅度
magnitude_sfm_noise_db =20*log10(magnitude_sfm_noise); %为画对数谱做准备
%以下解调
demodulated_signal_noise = fmdemod(s_fm_noise, fc, Fs_new, kfm); %调用自带解调函数，由于信号归一化，最大频偏就是kfm
demodulated_signal_resampled_noise = resample(demodulated_signal_noise, 1, 68);%重采样恢复原始音频信号
demodulated_signal_resampled_noise = demodulated_signal_resampled_noise / max(abs(demodulated_signal_resampled_noise)); %归一化
% audiowrite([file_path,'\测试解调效果_噪声.wav'], demodulated_signal_resampled_noise, Fs); %解调的声音存储为文件
%再计算解调后信号的频谱
y_de_noise =fft(demodulated_signal_resampled_noise); %对解调后的信号进行傅里叶变换
magnitude_de_noise = abs(y_de_noise(1:round(N/2))); %取幅度
magnitude_de_noise_db=20*log10(magnitude_de_noise); %为画对数谱做准备
%-------------------以下分析加噪声、使用预加重调制信号的频谱----------------
y_new_pre_noise = fft(s_fm_pre_noise); %fft
magnitude_sfm_pre_noise = abs(y_new_pre_noise(1:round(N_resampled/2))); %求幅度
f_new_pre = Fs_new/N_resampled*(0:round(N_resampled/2)-1); %显示实际频点的一半，频域映射，转化为HZ
magnitude_sfm_pre_noise_db =20*log10(magnitude_sfm_pre_noise); %为画对数谱做准备
%以下解调
demodulated_signal_pre_beforede_noise = fmdemod(s_fm_pre_noise, fc, Fs_new, kfm); %调用自带解调函数，由于信号归一化，最大频偏就是kfm
demodulated_signal_resampled_pre_beforede_noise = resample(demodulated_signal_pre_beforede_noise, 1, 68);%重采样恢复原始音频信号
demodulated_signal_resampled_pre_beforede_noise = demodulated_signal_resampled_pre_beforede_noise / max(abs(demodulated_signal_resampled_pre_beforede_noise)); %归一化
demodulated_signal_resampled_pre_noise = filter(1, [1 -0.9495], demodulated_signal_resampled_pre_beforede_noise);%去加重
% demodulated_signal_resampled_pre = resample(demodulated_signal_pre, 1, 68);%重采样恢复原始音频信号
demodulated_signal_resampled_pre_noise = demodulated_signal_resampled_pre_noise / max(abs(demodulated_signal_resampled_pre_noise)); %归一化
% audiowrite([file_path,'\预加重但未加重测试解调效果_噪声.wav'], demodulated_signal_resampled_pre_beforede_noise, Fs); %解调的声音存储为文件
% audiowrite([file_path,'\预加重去加重测试解调效果_噪声.wav'], demodulated_signal_resampled_pre_noise, Fs); %解调的声音存储为文件
%再计算解调后信号的频谱
% 去加重前
y_de_pre_beforede_noise =fft(demodulated_signal_resampled_pre_beforede_noise); %对解调后的信号进行傅里叶变换
magnitude_de_pre_beforede_noise = abs(y_de_pre_beforede_noise(1:round(N/2))); %取幅度
magnitude_de_pre_beforede_noise_db=20*log10(magnitude_de_pre_beforede_noise); %为画对数谱做准备
%去加重后
y_de_pre_noise =fft(demodulated_signal_resampled_pre_noise); %对解调后的信号进行傅里叶变换
magnitude_de_pre_noise = abs(y_de_pre_noise(1:round(N/2))); %取幅度
magnitude_de_pre_noise_db=20*log10(magnitude_de_pre_noise); %为画对数谱做准备
%==================比较用绘图=============
% ==================figure1===============================
figure;
subplot(3,1,1);%画三行一列的图
plot(t,x,'b');%绘制发送音频信号时域
hold on;
plot(t,x_filtered,'r');%绘制发送音频信号时域
hold off;
xlabel('Time/s');ylabel('Magnitude');
title('音频信号时域波形');
grid;
legend('滤波前','滤波后')
subplot(3,1,2);
plot(f,magnitude,'b'); %绘制发送音频信号频域
hold on;
plot(f,magnitude_filtered,'r'); %绘制发送音频信号频域
hold off;
xlabel('Frequency/Hz');ylabel('Magnitude');
title('音频信号频谱');
grid;
legend('滤波前','滤波后')
subplot(3,1,3);
plot(f,magnitude_db,'b'); %绘制原音频信号频域（对数谱）
hold on;
plot(f,magnitude_filtered_db,'r'); %绘制原音频信号频域（对数谱）
hold off;
xlabel('Frequency/Hz');ylabel('Magnitude/dB');
title('音频信号频谱（对数谱）');
grid;
legend('滤波前','滤波后')
%==================figure2===============================
figure;
subplot(3,1,1);%画三行一列的图
plot(t,x_filtered,'b');%绘制发送音频信号时域
hold on;
plot(t,x_preemph,'r');%绘制发送音频信号时域
hold off;
xlabel('Time/s');ylabel('Magnitude');
title('音频信号时域波形');
grid;
legend('预加重前','预加重后')
subplot(3,1,2);
plot(f,magnitude_filtered,'b'); %绘制发送音频信号频域
hold on;
plot(f,magnitude_preemph,'r'); %绘制发送音频信号频域
hold off;
xlabel('Frequency/Hz');ylabel('Magnitude');
title('音频信号频谱');
grid;
legend('预加重前','预加重后')
subplot(3,1,3);
plot(f,magnitude_filtered_db,'b'); %绘制原音频信号频域（对数谱）
hold on;
plot(f,magnitude_preemph_db,'r'); %绘制原音频信号频域（对数谱）
hold off;
xlabel('Frequency/Hz');ylabel('Magnitude/dB');
title('音频信号频谱（对数谱）');
grid;
legend('预加重前','预加重后')
%==================figure3===============================
figure;
subplot(3,1,1);%画三行一列的图
plot(t,demodulated_signal_resampled,'b');%绘制解调音频信号时域
hold on;
plot(t,demodulated_signal_resampled_pre,'r');%绘制解调音频信号时域
hold off;
xlabel('Time/s');ylabel('Magnitude');
title('未加噪声时解调信号时域波形');
grid;
legend('未使用预/去加重','使用')
subplot(3,1,2);
plot(f,magnitude_de,'b'); %绘制解调信号频域
hold on;
plot(f,magnitude_de_pre,'r'); %绘制解调信号频域
hold off;
xlabel('Frequency/Hz');ylabel('Magnitude');
title('未加噪声时解调信号频谱');
grid;
legend('未使用预/去加重','使用')
subplot(3,1,3);
plot(f,magnitude_de_db,'b'); %绘制解调信号频域（对数谱）
hold on;
plot(f,magnitude_de_pre_db,'r'); %绘制解调信号频域（对数谱）
hold off;
xlabel('Frequency/Hz');ylabel('Magnitude/dB');
title('未加噪声时解调信号频谱（对数谱）');
grid;
legend('未使用预/去加重','使用')
%==================figure4===============================
figure;
subplot(3,1,1);%画三行一列的图
plot(t_resampled,s_fm,'b');%绘制调制音频信号时域
hold on;
plot(t_resampled,s_fm_pre,'r');%绘制调制音频信号时域
hold off;
xlabel('Time/s');ylabel('Magnitude');
title('未加噪声时调制信号时域波形');
grid;
legend('未使用预/去加重','使用')
subplot(3,1,2);
plot(f_resampled,magnitude_sfm,'b'); %绘制调制信号频域
hold on;
plot(f_resampled,magnitude_sfm_pre,'r'); %绘制调制信号频域
hold off;
xlabel('Frequency/Hz');ylabel('Magnitude');
title('未加噪声时调制信号频谱');
grid;
legend('未使用预/去加重','使用')
subplot(3,1,3);
plot(f_resampled,magnitude_sfm_db,'b'); %绘制调制信号频域（对数谱）
hold on;
plot(f_resampled,magnitude_sfm_pre_db,'r'); %绘制调制信号频域（对数谱）
hold off;
xlabel('Frequency/Hz');ylabel('Magnitude/dB');
title('未加噪声时调制信号频谱（对数谱）');
grid;
legend('未使用预/去加重','使用')
%==================figure5===============================
figure;
subplot(3,1,1);%画三行一列的图
plot(t_resampled,s_fm_noise,'b');%绘制调制音频信号时域
hold on;
plot(t_resampled,s_fm_pre_noise,'r');%绘制调制音频信号时域
hold off;
xlabel('Time/s');ylabel('Magnitude');
title('dB噪声时调制信号时域波形');
grid;
legend('未使用预/去加重','使用')
subplot(3,1,2);
plot(f_resampled,magnitude_sfm,'b'); %绘制调制信号频域
hold on;
plot(f_resampled,magnitude_sfm_pre_noise,'r'); %绘制调制信号频域
hold off;
xlabel('Frequency/Hz');ylabel('Magnitude');
title('dB噪声时调制信号频谱');
grid;
legend('未使用预/去加重','使用')
subplot(3,1,3);
plot(f_resampled,magnitude_sfm_db,'b'); %绘制调制信号频域（对数谱）
hold on;
plot(f_resampled,magnitude_sfm_pre_noise_db,'r'); %绘制调制信号频域（对数谱）
hold off;
xlabel('Frequency/Hz');ylabel('Magnitude/dB');
title('dB噪声时调制信号频谱（对数谱）');
grid;
legend('未使用预/去加重','使用');
%==================figure6===============================
figure;
subplot(3,1,1);%画三行一列的图
plot(t,demodulated_signal_resampled_noise,'b');%绘制解调音频信号时域
hold on;
plot(t,demodulated_signal_resampled_pre_noise,'r');%绘制解调音频信号时域
hold off;
xlabel('Time/s');ylabel('Magnitude');
title('dB噪声时解调信号时域波形');
grid;
legend('未使用预/去加重','使用')
subplot(3,1,2);
plot(f,magnitude_de_noise,'b'); %绘制解调信号频域
hold on;
plot(f,magnitude_de_pre_noise,'r'); %绘制解调信号频域
hold off;
xlabel('Frequency/Hz');ylabel('Magnitude');
title('dB噪声时解调信号频谱');
grid;
legend('未使用预/去加重','使用')
subplot(3,1,3);
plot(f,magnitude_de_noise_db,'b'); %绘制解调信号频域（对数谱）
hold on;
plot(f,magnitude_de_pre_noise_db,'r'); %绘制解调信号频域（对数谱）
hold off;
xlabel('Frequency/Hz');ylabel('Magnitude/dB');
title('dB噪声时解调信号频谱（对数谱）');
grid;
legend('未使用预/去加重','使用')
%==================figure7===============================
figure;
subplot(3,1,1);%画三行一列的图
plot(t,demodulated_signal_resampled_pre,'b');%绘制解调音频信号时域
hold on;
plot(t,demodulated_signal_resampled_pre_beforede,'r');%绘制解调音频信号时域
hold off;
xlabel('Time/s');ylabel('Magnitude');
title('未加噪声时解调信号时域波形');
grid;
legend('使用预/去加重','使用预加重，未去加重')
subplot(3,1,2);
plot(f,magnitude_de_pre,'b'); %绘制解调信号频域
hold on;
plot(f,magnitude_de_pre_beforede,'r'); %绘制解调信号频域
hold off;
xlabel('Frequency/Hz');ylabel('Magnitude');
title('未加噪声时解调信号频谱');
grid;
legend('使用预/去加重','使用预加重，未去加重')
subplot(3,1,3);
plot(f,magnitude_de_pre_db,'b'); %绘制解调信号频域（对数谱）
hold on;
plot(f,magnitude_de_pre_beforede_db,'r'); %绘制解调信号频域（对数谱）
hold off;
xlabel('Frequency/Hz');ylabel('Magnitude/dB');
title('未加噪声时解调信号频谱（对数谱）');
grid;
legend('使用预/去加重','使用预加重，未去加重')
%==================figure8===============================
figure;
subplot(3,1,1);%画三行一列的图
plot(t,demodulated_signal_resampled_pre_noise,'b');%绘制解调音频信号时域
hold on;
plot(t,demodulated_signal_resampled_pre_beforede_noise,'r');%绘制解调音频信号时域
hold off;
xlabel('Time/s');ylabel('Magnitude');
title('信道信噪比dB时解调信号时域波形');
grid;
legend('使用预/去加重','使用预加重，未去加重')
subplot(3,1,2);
plot(f,magnitude_de_pre_noise,'b'); %绘制解调信号频域
hold on;
plot(f,magnitude_de_pre_beforede_noise,'r'); %绘制解调信号频域
hold off;
xlabel('Frequency/Hz');ylabel('Magnitude');
title('信道信噪比dB时解调信号频谱');
grid;
legend('使用预/去加重','使用预加重，未去加重')
subplot(3,1,3);
plot(f,magnitude_de_pre_noise_db,'b'); %绘制解调信号频域（对数谱）
hold on;
plot(f,magnitude_de_pre_beforede_noise_db,'r'); %绘制解调信号频域（对数谱）
hold off;
xlabel('Frequency/Hz');ylabel('Magnitude/dB');
title('信道信噪比dB时解调信号频谱（对数谱）');
grid;
legend('使用预/去加重','使用预加重，未去加重')
%==================figure9===============================
figure;
subplot(3,1,1);%画三行一列的图
plot(t_resampled,s_fm,'b');%绘制调制音频信号时域
hold on;
plot(t_resampled,s_fm_noise,'r');%绘制调制音频信号时域
hold off;
xlabel('Time/s');ylabel('Magnitude');
title('调制信号时域波形（未使用预加重）');
grid;
legend('无噪声','信噪比dB')
subplot(3,1,2);
plot(f_resampled,magnitude_sfm,'b'); %绘制调制信号频域
hold on;
plot(f_resampled,magnitude_sfm_noise,'r'); %绘制调制信号频域
hold off;
xlabel('Frequency/Hz');ylabel('Magnitude');
title('调制信号频谱（未使用预加重）');
grid;
legend('无噪声','信噪比dB')
subplot(3,1,3);
plot(f_resampled,magnitude_sfm_db,'b'); %绘制调制信号频域（对数谱）
hold on;
plot(f_resampled,magnitude_sfm_noise_db,'r'); %绘制调制信号频域（对数谱）
hold off;
xlabel('Frequency/Hz');ylabel('Magnitude/dB');
title('调制信号频谱（对数谱）（未使用预加重）');
grid;
legend('无噪声','信噪比dB')
%==================figure10===============================
figure;
subplot(3,1,1);%画三行一列的图
plot(t_resampled,s_fm_pre,'b');%绘制调制音频信号时域
hold on;
plot(t_resampled,s_fm_pre_noise,'r');%绘制调制音频信号时域
hold off;
xlabel('Time/s');ylabel('Magnitude');
title('调制信号时域波形（使用预加重）');
grid;
legend('无噪声','信噪比dB')
subplot(3,1,2);
plot(f_resampled,magnitude_sfm_pre,'b'); %绘制调制信号频域
hold on;
plot(f_resampled,magnitude_sfm_pre_noise,'r'); %绘制调制信号频域
hold off;
xlabel('Frequency/Hz');ylabel('Magnitude');
title('调制信号频谱（使用预加重）');
grid;
legend('无噪声','信噪比dB')
subplot(3,1,3);
plot(f_resampled,magnitude_sfm_pre_db,'b'); %绘制调制信号频域（对数谱）
hold on;
plot(f_resampled,magnitude_sfm_pre_noise_db,'r'); %绘制调制信号频域（对数谱）
hold off;
xlabel('Frequency/Hz');ylabel('Magnitude/dB');
title('调制信号频谱（对数谱）（使用预加重）');
grid;
legend('无噪声','信噪比dB')
关关难过关关过
u1s1
做了很久的心理建设 因为fengfeng之前的攻击力我是十分认可的
但是今天感觉还好
在我负责的部分
一个建议是预加重那个地方应该放原始音频
第二个是关于低信噪比时预加重效果反而差于普通组 温饱线以下高科技反而表现得更差 但是达到之后就表现得更好 人生哲理+1 让我再探究
第三个是关于 信噪比计算 解调增益计算
攻击力在承受范围之内
关于chirp信号过鉴频器 他也说挺好 It is my idea
汇报完之后心情还挺好的
标签：pre,noise,解调,When,Young,resampled,magnitude,信号,was From： https://www.cnblogs.com/LYoungH/p/17788225.html