clear;
clc;
close all;
%% 产生一个时域信号
SignalInput = cos(2 * pi * 1000 * (1:1:2048).' / 5000);
%% 设计一个带通滤波器
rp = 0.1; % Passband ripple
rs = 50; % Stopband ripple
fs = 2000; % Sampling frequency
f = [300 350 450 500]; % Cutoff frequencies
a = [0 1 0]; % Desired amplitudes
dev = [10^(-rs/20) (10^(rp/20)-1)/(10^(rp/20)+1) 10^(-rs/20)];
[n,fo,ao,w] = firpmord(f,a,dev,fs);
PassFreqcoeff = firpm(n,fo,ao,w); % 计算得到FIR带通滤波器的系数
%% 给信号添加带限高斯白噪声
Snr = 3;
[SignalAddNoise,Noisenor2] = SignalAddPassBandNoise(SignalInput,Snr,PassFreqcoeff);
%% 信噪比验证
ESignal = SignalInput' * SignalInput; % 信号能量
ENoise2 = Noisenor2 * Noisenor2'; % 噪声能量
aa = 10 * log10(ESignal/ENoise2); % 信噪比验证
%% 添加噪声后的信号
figure
plot(real(SignalAddNoise));
grid on
其中SignalAddPassBandNoise函数代码如下:
function [SignalAddNoise,Noisenor2] = SignalAddPassBandNoise(SignalInput,Snr,PassFreqcoeff)
%{
函数功能说明:给一段信号添加带限高斯白噪声;
输入参数说明:
SignalInput: 需要添加噪声的信号
Snr: 添加噪声后的信噪比
PassFreqcoeff: 带通滤波器系数
%}
%% 输入参数判断
% 输入的序列必须是N×1或者1×N;
[m,n] = size(SignalInput);
if (m~=1) && (n~=1)
error('输入信号SignalInput必须为向量形式!')
elseif m == 1
SendSignalNum = n;
TempSignalInput = SignalInput;
else
SendSignalNum = m;
TempSignalInput = SignalInput.';
end
%% 给信号添加噪声
Noise = randn(1, 2 * SendSignalNum); % 产生噪声点数是信号噪声的两倍
NoiseFilterd = filter(PassFreqcoeff, 1, Noise);
NoiseExtra = NoiseFilterd(SendSignalNum - SendSignalNum/2:SendSignalNum + SendSignalNum/2 - 1);
ENoise = NoiseExtra * NoiseExtra'; % 计算噪声的能量
ESignal = TempSignalInput * TempSignalInput'; % 计算信号的能量
Noisenorl = NoiseExtra / sqrt(ENoise); % 噪声归一化
Snr_ratio = 10^( - Snr/20);
Noisenor2 = Noisenorl * Snr_ratio * sqrt(ESignal); % 得到带限的高斯白噪声
SignalAddNoise = TempSignalInput + Noisenor2;
end
这个函数只能一维时域信号添加带限噪声,如果输入信号是一个矩阵则该函数无法实现。