任意频谱结构噪声生成

采用基于频率采样的FIR滤波器(MATLAB中为fir2)，生成任意频谱结构的噪声。

要点为 (1)频率的归一化; (2)功率谱密度的计算；(3)幅度的统一。

示例如下

% 设计任意形状的fir滤波器
% 有色噪声生成
clc; clear; close all
% 滤波器形状设计
ord = 48;  % 滤波器阶数
fs = 4000;
f_li = linspace(0, fs/2, ord);
mhi = zeros(1, ord);
mm = 0;
for i = 1: ord
    if f_li(i) <= 1000
        mhi(i) = f_li(i) * 5 + 50;
        temp = mhi(i);
    end 
    if f_li(i) > 1000
        mhi(i) = temp - mm;
        mm = mm + 10;
    end
end
mhi = mhi./100;

figure(1); 
subplot(211); plot(f_li, mhi); xlabel('频率/Hz'); ylabel('幅度'); grid on
subplot(212); plot(f_li, 20*log10(mhi)); xlabel('频率/Hz'); ylabel('幅度/dB'); grid on

f_li_1 = (f_li-min(f_li)) / (max(f_li)-min(f_li));  % 频率归一化到0~1
b = fir2(ord, f_li_1, mhi);                % 生成滤波器

figure(2); 
freqz(b, 1);
figure(2); hold on; 
subplot(211); plot(f_li/(fs/2), 20*log10(mhi), 'r');

% 生成白噪声
N = 10000;   % 采样点数
t = (0: 1: N-1)/fs;
x = sqrt(1)*randn(N, 1);% 生成白噪声，功率为1
sum(x.^2)/N   % 白噪声的功率

% 滤波
y_filted = filter(b, 1, x);
sum(y_filted.^2)/N   % 有色噪声的功率

figure(3)
subplot(211)
plot(t, x); xlabel('时间/s'); ylabel('幅度'); grid on; %ys = ylim;
title('原始噪声波形图');
subplot(212)
plot(t, y_filted); xlabel('时间/s'); ylabel('幅度'); grid on; %ylim(ys);
title('变换后噪声波形图');

% 频域
f = (0: 1: N/2)*fs/N;
X = fft(x);
X = abs(X(1: N/2+1)).^2 * 2/N/fs;   % 计算功率谱密度
Y_filted = fft(y_filted);
Y_filted = abs(Y_filted(1: N/2+1)).^2 * 2/N/fs;

figure(4)
subplot(211)
plot(f, 10*log10(X)); xlabel('频率/Hz'); ylabel('功率谱密度/dB·Hz^{-1}'); grid on
title('白噪声功率谱密度');
subplot(212)
plot(f, 10*log10(Y_filted)); xlabel('频率/Hz'); ylabel('功率谱密度/dB·Hz^{-1}'); grid on
hold on
plot(f_li, 20*log10(mhi)+mean(10*log10(X)), 'r');
legend('仿真频谱', '理想频谱', 'Location', 'best');
title('固定频谱结构的功率谱密度');