Andrew Hulva Digital Filtering 8/3/2015

Digital Filtering: Short Essay First things first, what is a digital filter? “A digital filter is a system that performs mathematical operations on a sampled, discrete-time signal to reduce or enhance certain aspects of that signal… this is in contrast to the other major type of electronic filter, the analog filter, which is an electronic circuit operating on continuous-time analog signals” (Wikipedia, 2015). A filter is characterized by its transfer or difference function. Generally, the transfer function has the form:

The order of the filter is the greater of N or M (Wikipedia, 2015). The above function does nothing, for the author, in the way of showing how to use a digital filter when approaching it mathematically; its generality is confusing. Luckily, there is an example that makes the visualization easier (Wikipedia, 2015). The process can be summed as such: a ratio of two polynomials that describes the filter in the z-domain are transformed to the discrete-time domain to yield a function where the next output sample (y(n)) can be described using some number of previous inputs, previous outputs, and current input. Starting with this function:

It is then expanded:

And divided by the highest order of z to make the corresponding filter casual (“a system where the output depends on past and current inputs but not future inputs”) (Wikipedia, 2015).

Rearranging the terms and preparing for the transform:

Andrew Hulva Digital Filtering 8/3/2015 Then taking the inverse z transform:

Finally, solving for y[n]:

Where y[n – p] are past outputs, x[n – p] are past inputs, x[n] is the current input, and y[n] is the next output (Wikipedia, 2015). This final function is what one can use to calculate the affected dataset (differential equations is quite useful after all). IIR vs. FIR In digital filtering, there are two types of filters that one hears of regularly: infinite impulse response and finite impulse response. What are they and what are the differences? An infinite impulse response filter (IIR) is a filter whose response to an impulse signal (usually Delta Dirac or Kronecker) does not become exactly zero past a certain point, or is infinite in its time length. This is due to feedback, or the next output value being derived from previous output, which results in the decay that never is exactly zero. In practice, this response eventually approaches zero and is rounded. The above example is of IIR type. Advantages of the IIR filter are that it requires less memory to calculate and is more efficient due to the feedback it receives from known quantities. Disadvantages of the IIR filter are that they are more susceptible to noise as a result of the rounding error (to make them have finite length) being compounded (Iowegian International Corporation, 2015). A finite response filter (FIR) is a filter whose response to an impulse is finite because its output past a certain time is indefinitely zero. The length of an Nth-order FIR filter is N + 1 samples or, equivalently, the number of coefficients. Advantages of the FIR filter are that is requires no feedback, so it does not suffer from compounded rounding error, and it can be easily designed to be phase coherent. The greatest disadvantage is that computation time is much longer than for a similar IIR filter.

Filtering Now use this knowledge to filter data in MATLAB.

Andrew Hulva Digital Filtering 8/3/2015 For the signal, “XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX” was loaded in MATLAB using the command: load(‘XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX’, ‘-mat’); % Alternatively, load the ‘testblast.mat’ attached with this document

Contained in this .mat file are a number of variables, those of importance for this lesson are ‘sampleRate’, ‘sensitivity’, ‘waveform’, ‘wOffset’, and ‘wScale.’ To convert the waveform from voltage to pressure: audio = (double(waveform) * wScale + wOffset) * 1000 / sensitivity;

Plot this audio signal: t = linspace(0, length(audio)/sampleRate, length(audio)); plot(t, audio)

Zoom in on region of interest and add labels: axis([.45 .75 -5 5]) xlabel('Time (s)') ylabel('Pressure (Pa)') title('Pressure Time Series'

The figure window should look like…

The signal of interest is the underlying low frequency information; specifically, its peak value and location are desired. The method by which this information will be found is by determining the location of a sign change in slope over some time duration, but the noise will disrupt this algorithm.

Andrew Hulva Digital Filtering 8/3/2015 To filter the data, a moving average filter will be built. This filter acts as a low-order low-pass filter about some constant frequency (the constant in the numerator of the transfer function). In this case, the constant component will be 0 Hz (a value of 1). To filter the data, three input arguments will be used with the ‘filter’ function in MATLAB. These are: “b” (numerator coefficients), “a” (denominator coefficients), and “x” (signal to be filtered and is already ‘audio’). Generate the necessary variables: windowSize = 100; b = (1/windowSize)*ones(1,windowSize); a = 1;

From the general form of the transfer function, this filter can be realized as: 𝐻(𝑧) =

𝐵(𝑧) 0.01 + 0.01𝑧 −1 + 0.01𝑧 −2 + ⋯ + 0.01𝑧 −100 = 𝐴(𝑧) 1

Notice how the denominator is unity, or 1, and does not feedback. This means the filter is of FIR type. Transforming this into the discrete time domain… 𝑦(𝑛) =

1 (𝑥(𝑛) + 𝑥(𝑛 − 1) + ⋯ + 𝑥(𝑛 − 99)) 100

Apply this filter to the noisy signal and view the output overlaid with the noisy signal: y = filter(b,a,audio); plot(t, y); hold plot(t, audio); axis([.45 .75 -5 5]) xlabel('Time (s)') ylabel('Pressure (Pa)') title('Pressure Time Series')

Andrew Hulva Digital Filtering 8/3/2015 The figure now should look like…

This looks similar to a 3-D image because of the delay inherent in the digital filter. The delay can be calculated with the formula: 𝑁−1 (𝑠𝑎𝑚𝑝𝑙𝑒𝑠) 2 Where, D = filter delay (samples) N = Order of filter (dB/octave) 𝐷=

To compensate for this delay, the first D samples will be removed and D zeros will be added to the end of the signal: delay = round((windowSize - 1) / 2); yy = y; yy(1:delay) = []; yy = [yy; zeros(delay,1)];

Andrew Hulva Digital Filtering 8/3/2015 Plot the two signals overlaid again: plot(t, audio) hold plot(t, yy) xlabel('Time (s)') ylabel('Pressure (Pa)') title('Pressure Time Series')

This should yield...

The signal of interest is now extracted and further processing can be done!

Andrew Hulva Digital Filtering 8/3/2015 Additional Information What about the frequency or impulse response of the filter? The frequency response of the filter can be viewed with the commands (requires DSP Toolbox): [filtresp, angfreq] = freqz(b,a); regfreq = angfreq * sampleRate / (2*pi); h2 = figure(); plot(regfreq, 20*log10(abs(filtresp))) axis([0 10000 -80 0]) xlabel('Frequency (Hz)') ylabel('Attenuation (dB)') title('Filter Frequency Response')

The impulse response can be viewed with the commands (requires DSP Toolbox): [impresp,impt] = impz(b,a); h3 = figure(); plot(impt, impresp); xlabel('Time (samples)') ylabel('y(n)') title('Filter Impulse Response') %Since the filter does not feedback, the impulse response is the constants across duration

A condensed version of the code (to prevent redundancy) is attached with this document. Running this will generate three figures (overlaid pressure-time series data, frequency response of filter, and impulse response of filter).

Andrew Hulva Digital Filtering 8/3/2015

Bibliography Iowegian International Corporation. (2015). IIR Basics FAQ. Retrieved from DSP Guru: http://dspguru.com/dsp/faqs/iir/basics Wikipedia. (2015, May 13). Casual System. Retrieved from Wikipedia, The Free Encyclopedia: https://en.wikipedia.org/wiki/Causal_system Wikipedia. (2015, August 3). Digital Filter. Retrieved from Wikipedia, The Free Encyclopedia: https://en.wikipedia.org/wiki/Digital_filter