4
$\begingroup$

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))));
Improve this question
$\endgroup$
7
  • $\begingroup$ How many values of the desired magnitude response do you have? It would be great if you could post (a link to) a figure showing the desired magnitude and your design result. $\endgroup$
    – Matt L.
    Jun 10, 2014 at 10:59
  • $\begingroup$ You can try to increase the order of your filter. A low order filter can fail to fit the desired frequency response. $\endgroup$
    – werediver
    Jun 10, 2014 at 11:02
  • $\begingroup$ @MattL. At the beginning I was using a polynomial to fit to my data and then design and filter accordingly. Then I also tried to use the exact signal. I will come back with an image. $\endgroup$
    – Christine Ioannidou
    Jun 10, 2014 at 11:21
  • $\begingroup$ @MattL.: dl.dropboxusercontent.com/u/43013925/filter%20response.png $\endgroup$
    – Christine Ioannidou
    Jun 10, 2014 at 12:00
  • $\begingroup$ OK, I see. It looks like the filter order is actually too high for the given number of data points. How many data points do you have? You could post the data (desired magnitude at specified frequencies), and I'll try to design a filter and show you how I would do it. $\endgroup$
    – Matt L.
    Jun 10, 2014 at 17:26

2 Answers 2

Reset to default
1
$\begingroup$

if you are only interested in designing a signal representing the wind turbine noise, you could just generate it from your magnitude. E.g., with your magnitude $A$ over frequency $f$ you can generate a random phase $\phi$ for each frequency and over a time $t$ and then add all together. Here a small code for illustration (I tried to use your varibles. However you need to check orientation etc.):

phi         = rand(size(fmag))*2*pi;
x           = (mag*sin(2*pi*fmag*t+repmat(phi,1,N)));

Hope this helps.

Improve this answer
$\endgroup$
0
$\begingroup$

Does remez() help you? Matlab has a similar one. This function receives some points of the spectrum response (PSD) and returns the coefficients of a filter with the desired response.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.