给信号添加指定信噪比的带限白噪声

1. 确定信号

首先我们生成确定信号。有了采样频率，我们需要设定采样点数，也是信号画图的点数N：N=2^n,信号点数一般取2的整数次方。

fo=800;       %中心频率800Hz
fs=40e3;      %采样频率40KHz
N=4096;       %采样点数
T=N/fs;       %信号长度,单位为s
t=(0:N-1)/fs; %信号时间序列
s1=sin(2*pi*fo*t);%正弦信号
figure(1)
subplot(221)
plot(t,s1);title('sin(2*pi*fo*t),fo=800Hz,fs=40KHz');
xlabel('时间/s');ylabel('幅度');%单频信号时域图

使用sin（）函数可以生成正弦信号，subplot（221）可以生成一个2*2的图阵，本图放在第一个。

线性调频信号是频率线性变化的信号，信号频率变化斜率为常数。

由于这是复信号模式，也可以只写它的实部。其中，K是信号频率变化斜率：

B=1e3;%带宽1KHz
k=B/T;%线性调频信号调频斜率
s2=exp(1j*(2*pi*fo*t+pi*k*t.^2));%线性调频信号从800Hz增加到1800Hz
subplot(223)
plot(t,real(s2));title('exp(j*(2*pi*fo*t+pi*k*t.^2)),fo=800Hz,fs=40KHz');
xlabel('时间/s');ylabel('幅度');%线性调频信号

对二者进行频谱分析，频谱分析使用FFT函数进行快速傅里叶变换。注意，FFT函数变换后，其幅值使原来的N/2倍。运行结果如下：

2. 随机信号
   为了便于观察信噪比，我设置随机信号长度为确定信号长度的3倍，使用normrnd函数生成指定均值方差的高斯随机序列：

N1=3*N;%噪声要较信号长，适合观察
mu=2;sigma=sqrt(2.5);%均值为2，方差为2.5
n=normrnd(mu,sigma,[1,N1]);%生成随机序列

带限白噪声是指在频域上限制其频带宽度。生成带限白噪声的方法有很多，可以使用时域加窗，比如海明窗、汉宁窗等，设计滤波器即可。这里我使用的是上课时学到的另一位同学的思路，他使用的是在频域直接加上矩形窗进行截断，这样生成的带限白噪声在频域十分简洁漂亮，便于观察。
我们直接将前面生成的白噪声转换到频域，对其加一个矩形窗，使其在500-2000Hz范围的幅值为1，其余为0.

fft_n1=[zeros(1,fix(N1/fs*fl)),fft_n(fix(N1/fs*fl):fix(N1/fs*fh)),zeros(1,fix(N1-N1/fs*fh))];%加矩形窗截断

其中，fix为取整函数，这是为了满足该向量赋值需要整数点。fft_n为我们转换到频域的噪声。fl,fh分别为噪声下限截止频率和上限截止频率。
这是运行结果，可以看出，对噪声进行带限处理后，明显噪声时域波形变得稀疏，这是因为其中500Hz以下和2000Hz以上的频率分量被滤掉了。

3. 规定谱级信噪比下的信号加噪声
   前面已经讲过信噪比和谱级信噪比的概念，所以这里我们直接进行公式推导：
   
   由此可以推出所需要的噪声能量：
   
   又因为采样后的信号是离散信号，使用sum（）求和代替积分，所以信噪比和噪声能量为：

E_n1=sum((real(n2)).^2);%噪声的能量
SNR1=10;
SNR2=0;
SNR3=-10;%三个谱级信噪比

求得所需要的噪声能量，又求得现在的噪声能量，对噪声幅值进行变化即可更改信号能量，得到所需信噪比：

s1_in=[zeros(1,N),s1(1:N),zeros(1,N)];%将信号与噪声统一长度
A_s1_SNR1=sqrt(P_s1/(10^(SNR1/10))*B_n);
n3_SNR1_s1=n2/sqrt(E_n1)*A_s1_SNR1;%与正弦信号达成谱级信噪比SNR1的噪声
s1_SNR1_out=s1_in+n3_SNR1_s1;%正弦信号加入噪声后的输出信号

运行结果为：

如果内容有什么问题，欢迎大家指导我，matlab新手，状况百出hhhh