1/3oct smoothing 是什么意思，怎样实现 - 基础理论室 - 声学楼论坛

下段供参考吧,注意中间的smooth的一段.

我想,你想一想就能理解的.

我的理解是在那段频带内(1/3或1/6倍频程)的数据求平均,

具体的算法要查那些仪器和测试系统的相关资料.

原则上就是你已经查到的那些算法中的某一种.

过去手算的时侯,就是求算术平均值.

-----------

Convolution

In the real world it would be nice to do Fourier transform filtering; however, it takes a lot of computing power to do a FFT. If one is working offline (e.g. turning and audio file into an mpeg encoded file) the time taken for Fourier processing does not really matter. For signals that must be processed in real time (e.g. live audio) the computing overhead required to calculate an FFT do the filtering operation and then convert back with an IFFT becomes impractical either because the process cannot be done fast enough or because the computing power is too expensive. It turns out that there is a way of designing a filter in the Fourier domain and then implementing that filter digitally in the time domain. The process of using the filters time domain version involves a mathematical process called convolution. I’ll explain convolution in class, then we will do this little exercise to appreciate what a convolution does to simple wave forms.

1. To see what convolution does enter a square pulse signal, x, consisting of 10 0’s 10 1’s and 10 0’s. Like this.

>> x=[0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0]

Plot this signal. What do you expect when you convolve this signal with itself? You should be able to figure out what the result will look like. Now try to form the convolution in practice using MATLAB’s conv function.

>> xx=conv(x,x);

Plot xx the result of convolving x with itself. Is it what you predicted?

Now try a different function—the triangle

>> p=[0 0 0 0 1 2 3 4 5 4 3 2 1 0 0 0 0]

>>pp=conv(p,p);

Plot it. Pretty, no?

2. It is not necessary to convolve a function only with itself. In the next example we will use convolution to reduce the noise in a signal. First let’s create a noisy signal. Create a sine wave

>> t=0:0.001:1;

>> w=2*pi*3;

>> y=sin(w*t);

Plot it. It is nice and “clean”. What is the sampling rate in my example? Does the plot make sense in terms of number of waves in the 1 second time window?

Now we create noise using the rand function.

>> for k=1:max(size(t))

noise(k)=rand*0.2;

end

Next we add the noise to the original signal

>>yn=y+noise;

Plot the array of signal plus noise yn. Hold that plot!

The smoothing function that we will use is based on binomial coefficients obtained from Pascal’s triangle. Again some class explanation is necessary to understand how to create these coefficients and (more importantly) why they are useful for smoothing.

Create the binomial smoothing function

>> binom=(1/64)*[1 6 15 20 15 6 1]

Why do we multiply by (1/64)?

Smooth the signal by convolving the noisy signal yn with the smoothing function.

yfilt=conv(binom,yn);

Plot the filtered yfilt on top of the signal plus noise yn. Smoother, no? Try applying the filter again. Plot. Even more smooth, no? [This text sound like it was written by Pepe le Pew!] Why does the plot seem to shift to the right with each smoothing?

3. Now for the real deal—creating a filter in the time domain by specifying the filter in the frequency domain and doing an inverse Fourier transform. Our aim is to create a low pass filter that has a cutoff of 500 Hz. We will test this filter by creating complex tone with a bunch of harmonics some of which go beyond the pass band.

Step 1—creating the filter in the frequency domain. Our sampling rate will be our old favorite 22050. In the frequency domain our ideal low pass filter has the first 500 Hz equal to 1 and all the rest equal to 0. We call the frequency domain filter profile H. This profile is referred to as the transfer function of the filter.

>> H=zeros(1,22051);

>> H(1:500)=1;

>> H(22051-500:22051)=1;

Step2—turn this profile in frequency into a time domain filter function. The time domain signal is called the impulse response.

>> filt=ifft(H);

Now convert this impulse response so that it is correctly displayed in time and then plot it. OK I know some explanation is needed on the fftshift idea…

>>filt=fftshift(filt);

>>plot(real(filt))

Look carefully at the impulse response. First, it goes on across all time values. In reality it goes on in time from negative infinity to positive infinity. Such filter behavior is called an Infinite Impulse Response (IIR). Clearly such infinite behavior will be a bitch to deal with in practice. Luckily the signal drops off in each direction so that we can truncate the filter and only use a small section. Although this point is not clear from our analysis, the impulse goes on from times that precede t=0. Such behavior is called acausal. It implies that filtering starts before any signal is present. To create our truncated signal, h, we want 500 points on either side of the central maximum.

>>h=real(filt(10526:11526));

Step 3—lets test our filter. Create a complex tone of fundamental 200 Hz and a bunch of harmonics (I chose odd harmonics for no particularly good reason.

>> w=2*pi*200

>> y=sin(w*t)+sin(3*w*t)+sin(5*w*t)+sin(7*w*t)+sin(9*w*t)+sin(11*w*t);

Form the spectrum of our complex tone and plot it and hold that plot.

>> fy=fft(y);

>> plot(abs(fy(1:11025))

>>hold

Now we use our time domain representation of the filter, h, and apply it to our complex tone, y, to get our filtered signal, fsig.

>> fsig=conv(h,y);

Take the spectrum of fsig and plot it on top of the spectrum of our original signal (in red). Note that when I take the FFT I cut off the fist 500 and the last 500 elements of fsig. Why? How come fsig is longer than y?

>> ffsig=fft(fsig(500:22550));

>> plot(abs(ffsig(1:11025)),'r')

Has filtering taken place in the manner we expected? Did this process remove the higher harmonics? Which ones did it remove?

4. Final convolution exercise. The last filtering exercise was performed on a signal with a bunch of well spaced apart frequencies. That example made the filter look pretty good. To appreciate that there are some limitations we now use our same filter from last time, h, applied to a continuous set of frequencies. To get a continuous frequency distribution we use our old friend click.m to generate a short “click”. As we have seen in the past the click possesses a broad continuous range of frequencies. Choose s=0.0001 (I think? We might experiment in class. Run click. Plot the click, y, to remind yourself of the time signal we are using. Take an FFT of y and plot and hold it.

>> click

>> plot(y)

>> fy=fft(y);

>> plot(abs(fy(1:11025)));

>> hold

Now filter the click using our impulse response filter, h.

>>yfilt=conv(h,y);

Form the FFT of our filtered signal yfilt and plot it.

>> fyfilt=fft(yfilt(500:22550));

>> plot(abs(fy(1:11025)))

Zoom in to see if the cut-off is nice and clean. Is it perfect? Why not?

5. Windowing. Using the truncated time signal, h, (only 1000 time steps) instead of the full time signal, filt, (22050 time steps) leads to a frequency cut-off that has a bunch of little bumps beyond our desired cut-off frequency. This problem is very common in time and spectrum manipulation. It occurs whenever the signal in either time or frequency is cut off suddenly. Nature does not like that sharp, abrupt cut-off. For example, it is really not possible to synthesize a square wave that is perfectly square because it would require an infinite set of harmonic frequencies. If you think back to the square wave exercise we did earlier we always had a square wave with wiggly shoulders because we only had a finite number of harmonics in our synthesis. In the case here we truncated the time signal, h, by just cutting it off—essentially we multiplied by a square time window 1000 time steps long. Nature doesn’t like those square windows, and, as a consequence, we ended up with the bumps. We can get rid of those ugly bumps by the use of a smooth window instead of the abrupt square window. Let’s make the smooth window—there are many and this one is called the Blackman window. Open a new m-file and type the following

function [w]=blackmanwind(M)

% Generate a Blackmann window length M

w = .42 - .5*cos(2*pi*(0:M-1)/(M-1)) + .08*cos(4*pi*(0:M-1)/(M-1));

Save this file as blackmanwind.m

Back in the Command window run

>> [w]=blackmanwind(1001);

This creates w the new, non-square window. Plot it. Not square, eh!

Now multiply h by the Blackman window and plot it. See how the signal h2 tapers away smoothly to zero at the edges.

>> h2=h.*w;

>> plot(h2)

Now convolve our new time domain filter function h2 with the click signal y; form a spectrum and plot it.

>> convy2=conv(h2,y);

>> fconvy=fft(convy(500:22550);

Plot the bumpy signal from before and hold it.

>> plot(fyfilt);

>>hold

Now plot our Blackmann window spectrum in red.

>>plot(fconvy,’r’)

Examine the red curve compared to the blue—smoother, no?

Epilogue

Why are these time domain representations easier to use? Because they can be implemented on a signal as it passes through a digital signal processing (DSP) module. DSP’s have registers and buffers to hold a section of the signal and perform the time convolution as the signal passes through the device. This process is much easier and cheaper to implement compared to FFTs. Of course we just looked at one simple filter—the possibilities are enormous.

主题：1/3oct smoothing 是什么意思，怎样实现