When I was Young

%以下代码功能：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 

汇报完之后心情还挺好的