信息处理 - 如何设计具有一定幅度响应的滤波器 - 吾爱随笔录

我正在尝试设计一个幅度与给定信号相同的滤波器。给定的信号是风力涡轮机噪声，因此它具有显着的低频成分。设计滤波器后，我想过滤高斯白噪声，以创建风力涡轮机噪声模型。这两个信号，即原始噪声和过滤后的噪声应该听起来相似。

我在 Matlab 中使用任意幅度滤波器设计（FIR，阶数：900，单速率，1 波段，由幅度指定的响应，采样率 44100 Hz，设计方法：firls）。问题是，尽管我使用来自原始信号幅度的值来设计滤波器，但滤波器幅度未能遵循更高频率的幅度。你能帮我解决这个问题吗？

这个想法是我使用多项式曲线来拟合原始声音的频谱形状。然后，我通过将提取的多项式曲线与生成的噪声幅度相乘来过滤生成的噪声的幅度。最后，在计算出新的幅度和相位之后，我回到了时域。

[x,fs] = audioread('cotton_0115-0145.wav');  % original noise sample
x = x(:,1);                                  % extract one channel
x = x.';
N = length (x);

% fft of original signal
Y = fft(x,N)/N; 
fy = (0:N-1)'*fs/N;

% half-bandwidth selection for original signal
mag = abs(Y(1:N/2));
fmag = fy(1:N/2);

% polynomial fitting
degreesOfFreedom=40;  
tempMag=20*log10(mag)';  % desired magnitude in dB
tempFmag=fmag;

figure(1)
plot(tempFmag,tempMag,'b');
title('Spectrum of original signal-Polynomial fitting')
xlabel('Frequency (Hz)'); 
ylabel('20log10(abs(fft(x)))');
axis([tempFmag(1) tempFmag(end) min(tempMag) 0]);

hold on
p = [];

for i=1:4
    p=polyfit(tempFmag,tempMag,degreesOfFreedom);
    pmag=polyval(p,tempFmag);

    plot(tempFmag,pmag,'r');
    pause

    above=pmag<tempMag;
    abovemag=tempMag(above);
    abovefmag=tempFmag(above);

    tempMag=abovemag;
    tempFmag=abovefmag;

end

hold on
legend('signal magnitude','polynomial');

%
loc1 = find(fmag == 0);
loc2 = find(fmag == 22050);
Nmag = length(mag);

M=((Nmag-1)*max(tempFmag))/(tempFmag(end)-tempFmag(1));
freqFinal=tempFmag(1):max(tempFmag)/M:max(tempFmag);
freqFinal=tempFmag(1):max(tempFmag)/length(mag):max(tempFmag);
magnitudesFinal=polyval(p,freqFinal);

figure(2) 
plot(fmag,20*log10(mag)');
hold on;
plot(freqFinal,magnitudesFinal,'g');
title('Spectrum of original signal-Choice of polynomial curve')
xlabel('Frequency (Hz)'); 
ylabel('abs(fft(x))');
axis([freqFinal(1) freqFinal(end) min(magnitudesFinal) max(magnitudesFinal)]);

%%
% noise generation
Nn  = N;
noise =  wgn(1,Nn,0);
noise = noise/(max(abs(noise)));
Ynoise = fft(noise,Nn)/Nn;
fn = (0:Nn-1)'*fs/Nn;

% polynomial for whole f range
newmagA = 10.^(magnitudesFinal/20);
newmagB = fliplr(newmagA);
newmagB(end+1) = newmagB(end);
newmag = [newmagA newmagB];

%filtering
Ynoisenew = newmag .* Ynoise;
figure(3)
magnoise = 20*log10(abs(Ynoisenew(1:Nn/2)));
fnoise = fn(1:Nn/2);
plot(fnoise, magnoise);

% magnitude and phase of filtered noise
magn = abs(Ynoisenew);
phn = unwrap(angle(Ynoisenew));

% Back to Time domain
sig_new=real(ifft(magn.*exp(1i*phn)));

figure(4)
sig_new = sig_new / max(abs(sig_new));
plot(t, sig_new);

Ysignew = fft(sig_new,Nn)/Nn;
fn = (0:Nn-1)'*fs/Nn;
figure(5); plot(fn(1:Nn/2),20*log10(abs(Ysignew(1:Nn/2))));