I am trying to design a filter whose magnitude is the same as that of a given signal. The given signal is wind turbine noise, so it has significant low-frequency content. After designing the filter, I want to filter white Gaussian noise so as to create a model of wind turbine noise. The two signals, that is the original and the filtered noise should sound similar.
I am using arbitrary magnitude filter design in Matlab for that (FIR, Order: 900, Single rate, 1-band, response specified by amplitudes, Sample rate 44100 Hz, Design method: firls). The problem is that, although I design the filter using the values from the original signal's magnitude, the filter magnitude fails to follow the magnitude at higher frequencies. Could you please help me with that?
The idea is that I am using a polynomial curve, to fit the shape of the spectrum of the original sound. Then, I filter the magnitude of the generated noise, by multiplying the extracted polynomial curve with the generated noise magnitude. Finally, after calculating the new magnitude and phase, I get back to the time domain.
[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))));