Sound Processing and Digital Filters Graduate School of Culture - - PowerPoint PPT Presentation

1

Sound Processing and Digital Filters

CTP 431 Music and Audio Computing

Graduate School of Culture Technology (GSCT) Juhan Nam

Outlines

§ Introduction: Sound Processing § Linear Time-Invariant Digital Filter

– Impulse response – Convolution

§ Digital Filters

– FIR Filters – IIR Filters

§ Frequency Response § Transfer functions

– Z-transform – Pole-Zero Analysis

§ Bi-quad Filters

2

Introduction: Sound Processing

§ Sounds captured on computers are processed in various ways

– Sample editing: cut, copy, paste – Amplitude: gain, fade in/out, automation curve, compressor – Timbre: lowpass/highpass filters, EQ, distortion, modulation – Spatial effect: delay, reverberation – Pitch: pitch shifting (e.g. auto-tune) – Time stretching – Noise suppression

3

Sound Processing

§ Linear processing

– No sinusoidal components are introduced by the processing

• Only the amplitude and phase of sinusoidal components of the

input change – Filters (lowpass, highpass, bandpass, …), EQ – Delay-based audio effect: delay, chorus, flanger, reverberation

§ Non-linear processing

– New sinusoidal components are generated – Compressor, distortion, clipping – Pitch shifting, ring modulation, …

4

Processor Input Output

Linear Time-Invariant (LTI) Digital Filters

§ What linear filters can do

– Amplifying the input (or the past output): e.g. y(n) = ax(n) – Delaying the input (or the past output): e.g. y(n) = x(n-1) – Summing them all: y(n) = ax(n) + x(n-1)

§ “Easy-to-understand” definition

5

Input Output

x(n) y(n)

Filter

LTI Digital Filters (Formal)

§ Linearity

– Homogeneity: if x(n) à y(n), ax(n) à ay(n) – Superposition: if x1(n) à y1 (n) and x2(n) à y2 (n), x1(n) + x2 (n) à y1(n) + y2 (n)

§ Time-Invariance

– If x(n) à y(n), x(n-N) à y(n-N) for any N – The system does not change its behavior with time.

• In practice, most systems do change over time but not quickly

6

Input Output

x(n) y(n)

Filter

Example: Simple LTI Digital Filters

§ Moving-average filter

– y(n) = 0.5x(n) + 0.5x(n-1) – Low-pass

§ Differentiator

– y(n) = 0.5x(n) - 0.5x(n-1) – High-pass

§ Feed-forward comb filter

– y(n) = x(n) + x(n-M) where M is, say,100 – Renders harmonically distributed peaks and valleys

7

Impulse Response

§ Characterize filters as a number sequence § Obtained when x(n) is a unit impulse

– x(n) = δ(n) = [1, 0, 0, 0, …] à y(n) = h(n)

§ Can be measured from a linear system (black-box approach)

– If you excite the linear system with an impulse, you can record the

• utput and use that to determine exactly what the system response

would be to any arbitrary input.

8

Input Output Filter

h(n)

x(n) = δ(n) y(n) = h(n)

Convolution

§ The output of LTI systems is represented by convolution

• peration between the input x(n) and impulse response

h(n) § Deriving convolution

– The input is decomposed into a time-ordered set of weighted impulse

• x(n) = [x0, x1, x2, x3, … , xM, ]

= x0δ(n) + x1δ(n-1) + x2δ(n-2) + x3δ(n-3) + …+ xMδ(n-M) – By the linearity and time-invariance

• y(n) = x0h(n) + x1h(n-1) + x2h(n-2) + x3h(n-3) + …+ xMh(n-M)

9

y(n) = x(n)*h(n) = x(i)h(n −i)

i=0 M

∑

Convolution in Practice

§ From the commutative law

– y(n) = h0x(n) + h1x(n-1) + h2x(n-2) + h3x(n-3) + …+ hMx(n-M) – For example: x(n) (n=0,1,2,3,4,5) , h(n) = [h0 h1 h2 ]

• y(0) = h0x(0)
• y(1) = h0x(1) + h1x(0)
• y(2) = h0x(2) + h1x(1) + h2x(0)
• y(3) = h0x(3) + h1x(2) + h2x (1)
• y(4) = h0x(4) + h1x(3) + h2x (2)
• y(5) = h0x(5) + h1x(4) + h2x (3)
• y(6) = h1x(5) + h2x (4)
• y(7) = h2x (5)

10

y(n) = x(n)*h(n) = x(i)h(n −i)

i=0 M

∑

= x(n −i)h(i)

i=0 M

∑

Fully overlapped region (steady-state) Transient region Transient region

Convolution in Practice

§ If the length of x(n) is M and the length of h(n) is N, the length of y(n) is M+N-1 § Computation Complexity

– In Big-O notation, it requires M x N multiplications – If N is a large number, it is quite expensive to compute – We can compute convolution in frequency domain, which is much cheaper than the time-domain approach when the impulse response is long (e.g. reverberation or head-related transfer functions)

11

Properties of LTI systems

§ Commutative § Associative § Distributive

12

x(n)*h

1(n)*h2(n) = x(n)*h2(n)*h 1(n)

{x(n)*h

1(n)}*h2(n) = x(n)*{h 1(n)*h2(n)}

x(n)*h

1(n)+ x(n)*h2(n) = x(n)*{h 1(n)+ h2(n)}

FIR Filters

§ The output is formed from input and its past input

– They have finite impulse responses (FIR) – Convolution with the finite impulse responses

13

“Delay operator” (a unit sample of delay)

Z-1 Z-1 x(n)

+

y(n) Z-1 Z-1 . . . b1 b2 bM b0 x(n-1) x(n-2) x(n-M)

y(n) = b0x(n)+ b

1x(n −1)+ b2x(n − 2)+...+ bM x(n − M)

IIR Filters

§ The output can be also formed by input and past outputs

– The feedback creates an infinite impulse response! – Convolution with the infinite impulse responses

14

Z-1 x(n)

+

y(n) Z-1 Z-1 . . .

• a1
• a2
• aN

b0 y(n-1) y(n-2) y(n-M)

y(n) = b0x(n)− a1y(n −1)− a2y(n − 2)−...− bNy(n − N)

IIR Filters

§ The infinite impulse response

– For example: y(n) = x(n) + ry(n-1)

• y(0) = x(0)
• y(1) = x(1) + ry(0) = x(1) + rx(0)
• y(2) = x(2) + ry(1) = x(2) + rx(1) + r2x(0)
• y(3) = x(3) + ry(2) = x(3) + rx(2) + r2x(1) + r3x(0)
• y(4) = x(4) + ry(3) = x(4) + rx(3) + r2x(2) + r3x(1) + r4x(0)

– Stability issue!

• If r < 1, the filter becomes stable
• If r =1, the filter oscillates
• If r > 1, the filter becomes unstable

– The impulse response is long but, in practice, it is finite (for r < 1)
 because the level goes below the the quantization noise floor

15

à h(n) = [1 r r2 r3 r4 r5 r6 …]

Example: IIR Filters

§ Leaky Integrator

– y(n) = x(n) + ry(n-1) – r is a slightly less than 1. (1-r) is the “leak” – Lowpass filtering

§ “Reson” filter

– y(n) = x(n) +2rcosθy(n-1) - r2y(n-2) – r controls the resonance and θ controls its frequency

• Resonance: 0 (low resonance) < r < 1 (low resonance)
• Cut-off frequency (fc): θ=2πfc/fs (fs: sampling rate)

– Low-pass/band-pass/high-pass filtering depending on additional zeros

§ Feed-back comb filter

– y(n) = x(n) + ry(n-M) – Renders a harmonic tone

16

General Filter Form

§ The general form of digital Filters

– The order of the filter is the greater of M or N

17

y(n) = b0x(n)+ b

1x(n −1)+ b2x(n − 2)+...+ bM x(n − M)

− a1y(n −1)− a2y(n − 2)−...− aNy(n − N)

Z-1 x(n)

+

Z-1 Z-1 . . . b1 b2 bM b0 x(n-1) x(n-2) x(n-M) Z-1 y(n) Z-1 Z-1 . . .

• a1
• a2
• aN

y(n-1) y(n-2) y(n-N) “Difference EquaLon”

Filter Implementation Forms

§ Direct Form I § Direct Form II

18

Z-1 x(n)

+

Z-1 b1 b0 Z-1 y(n) Z-1

• a1

b2

• a2

x(n)

+

Z-1 Z-1

• a1
• a2

+

b1 b0 b2 y(n) ( By the commutative rule) “Biquad Filter”

Example of Filter Implementation

§ Typically implemented in time-domain

19

A short version using “filter” function in Matlab

Frequency Responses

§ Describe the characteristics of filters in frequency domain

– Amplitude response: the amount of amplitude change (often in dB) – Phase response: the amount of delay (-2pi ~ 0)

20

10

2

10

3

10

4

−30 −20 −10 10 20 30 Amplitude Response freqeuncy(log10) Gain(dB) 10

2

10

3

10

4

−3 −2.5 −2 −1.5 −1 −0.5 Phase Response freqeuncy(log10) radian

Frequency Responses

21

Input Output Filter x(n) = ejωn y(n) = G(ω)ej(ω+Θ(ω))

20 40 60 80 100 120 140 160 180 200 −1 −0.5 0.5 1 Time [Samples] Amplitude Input Frequency = 100Hz input

• utput

20 40 60 80 100 120 140 160 180 200 −1 −0.5 0.5 1 Time [Samples] Amplitude Input Frequency = 1000Hz input

• utput

§ For the sinusoidal input and outputs

– x(n) = ejωn à x(n-m) = ejω(n-m) = e-jωm x(n) for any m – y(n) = G(ω)ej(ωn+Θ(ω)) à y(n-m) = G(ω)ej(ω(n-m)+Θ(ω)) = e-jωm y(n) for any m

§ Putting this property into the general form of difference equation

22

Frequency Response

y(n) = b0x(n)+ b

1e− jωx(n)+ b2e− j2ωx(n)+...+ bMe− jMωx(n)

− a1e− jωy(n)− a2e− j2ωy(n)−...− aNe− jNωy(n)

y(n) = b0 + b

1e− jω + b2e− j2ω +...+ bMe− jMω

1+ a1e− jω + a2e− j2ω +...+ aNe− jNω x(n) H(ω) = b0 + b

1e− jω + b2e− j2ω +...+ bMe− jMω

1+ a1e− jω + a2e− j2ω +...+ aNe− jNω

H(ω) : frequency response |H(ω)| = G(ω) : amplitude response H(ω) = Θ(ω) :phase response

∠

Examples of Frequency Response

§ Moving-average filter (lowpass)

– y(n) = 0.5(x(n) + x(n-1)) – H(ω) = 0.5(1+ e-jω) = 0.5(ejω/2+ e-jω/2)e-jω/2 = cos(ω/2) e-jω/2 – G(ω) = cos(ω/2), Θ(ω) = -ω/2

§ Differentiator (highpass)

– y(n) = 0.5(x(n) + x(n-1)) – H(ω) = 0.5(1- e-jω) = 0.5(ejω/2- e-jω/2)e-jω/2 = sin(ω/2) e-jω/2+jπ/2 – G(ω) = sin(ω/2), Θ(ω) = -ω/2+π/2

23

Z-Transform

§ Z-transform

– Define z to be a variable in complex plane: we call it z-plane – When z = ejω (on unit circle), the frequency response is a particular case of the following – We call this Z-transform of h(n) or transfer function – z-1 corresponds to one sample delay:

• Called delay operator or delay element

– Filters are often expressed as Z-transform

24

H(z) = B(z) A(z) = b0 + b

1z−1 + b2z−2 +...+ bMz−M

1+ a1z−1 + a2z−2 +...+ aNz−N

Poles and Zeros

§ H(z) can be factorized and we can find roots for each of polynomials

– Zeros: the numerator roots – Poles: the denominator roots

§ We can analyze frequency response of filters more easily with poles and zeros than numerator or denominator coefficient!!!

25

H(z) = B(z) A(z) = (1− q1z−1)(1− q2z−1)(1− q3z−1)...(1− qMz−1) (1− p1z−1)(1− p2z−1)(1− p3z−1)...(1− pNz−1)

Pole-Zero Analysis: Amplitude Response

§ The amplitude response is represented as

– Numerator: factors of distance between zero and unit circle – Denominator: factors of distance between pole and unit circle

26

G(ω) = H(z = e jω ) = (1− q1e− jω )(1− q2e− jω )(1− q3e− jω )...(1− qMe− jω ) (1− p1e− jω )(1− p2e− jω )(1− p3e− jω )...(1− pNe− jω ) = (e jω − q1)(e jω − q2)(e jω − q3)...(e jω − qM ) (e jω − p1)(e jω − p2)(e jω − p3)...(e jω − pN ) = (e jω − q1) (e jω − q2) (e jω − q3)... (e jω − qM ) (e jω − p1) (e jω − p2) (e jω − p3)... (e jω − pN )

Pole-Zero Analysis: Amplitude Response

§ Bi-quad case

– The amplitude response is given as – As poles are close to the unit circle,
 the amplitude response is boosted – As zeros are close to the unit circles,
 the amplitude response is damped

§ If poles are outside the unit circle, the filter becomes unstable!

– If poles are on the unit circles, the filter oscillates.

27

G(ω) = d1(ω)d2(ω) d3(ω)d4(ω)

Examples

§ Moving-average filter (lowpass)

– y(n) = 0.5(x(n) + x(n-1)) – zeros: z = -1 (no poles)

§ Leaky Integrator

– y(n) = x(n) + ry(n-1) – poles: z = -r (no zeros)

§ Reson filter

– y(n) = x(n) +2rcosθy(n-1) - r2y(n-2) – poles: z = r(cosθ + jsinθ), r(cosθ - jsinθ) (no zeros)

§ Feed-back comb filter

– y(n) = x(n) - ry(n-M) (for convenience, the sign of r has changed) – poles: z = r1/M (cos(2π/Mn) + jsin(2π/Mn)) (n=0, 1, 2, … M-1) (no zeros)

28

Pole-Zero Analysis: Phase Response

§ The phase response is represented as

29

Θ(ω) = ∠H(z = e jω ) = ∠ (1− q1e− jω )(1− q2e− jω )(1− q3e− jω )...(1− qMe− jω ) (1− p1e− jω )(1− p2e− jω )(1− p3e− jω )...(1− pNe− jω ) = ∠(e jω − q1)+∠(e jω − q2)+∠(e jω − q3)...∠(e jω − qM ) −∠(e jω − p1)−∠(e jω − p2)−∠(e jω − p3)...−∠(e jω − pN )

Pole-Zero Analysis: Phase Response

§ In the following examples, the phase response is given as

– Positive angles for zeros – Negative angle for poles

30

Θ(ω) =θ1 +θ2 −θ3 −θ4

Formal Definition of Z-transform

§ Z-transform § Properties

– Shift theorem: – Convolution theorem:

• Therefore, the transfer function is represented as

§ Decomposing z-transforms

– Series combination: – Parallel combination:

31

X(z) = x(n) n=0 ∞

∑

z−n x(n − Δ)↔ z−ΔX(z) x(n)*h(n)↔ X(z)H(z) y(n) = x(n)*h(n)↔ H(z) = Y(z) X(z) H(z) = H1(z)H2(z)↔ h(n) = h

1(n)*h2(n)

H(z) = H1(z)+ H2(z)↔ h(n) = h

1(n)+ h2(n)

Frequency Response by DTFT

§ Discrete-Time Fourier Transform

– Putting back to the Z-transform

§ The properties work for DTFT in the same manner

– Shift theorem: – Convolution theorem:

§ The frequency response is represented as

32

z−1 = e− jω X(e− jω ) = X(ω) = x(n) n=0 ∞

∑

e− jnω y(n) = x(n)*h(n)↔ H(ω) = Y(ω) H(ω) x(n − Δ)↔ e− jΔωX(ω) x(n)*h(n)↔ X(ω)H(ω)

Practical Filters

§ One-pole one-zero filters

– Moving average filters: low-pass filter – Leaky integrator: low-pass filter – DC-removal filters: high-pass filter – Bass / treble shelving filter

§ Bi-quad filters

– Low-pass / high-pass filters – Band-pass / botch filters – EQ – The two-poles are real-numbers or complex conjugate

§ Note that any high-order filter can be factored into one- pole one-zero filters and bi-quad filters!

33

H(z) = b0 + b

1z−1

a0 + a1z−1

H(z) = b0 + b

1z−1 + b2z−2

a0 + a1z−1 + a2z−2

One-pole one-zero filters

§ DC-removal filters: high-pass filter

34

H(z) = 1− z−1 1− az−1 (1+ a 2 )

10

2

10

3

10

4

−40 −35 −30 −25 −20 −15 −10 −5 5 10 a =0.9 a =0.92 a =0.94 a =0.96 a =0.98 DC removal filter: Amplitude Response freqeuncy(log10) Gain(dB)

One-pole one-zero filters

§ Bass / Treble shelving

35

H(z) =1+ g 1+ z−1 1− az−1 1− a 2 H(z) =1+ g 1− z−1 1− az−1 1+ a 2

g =10

GdB 20 −1 10

2

10

3

10

4

−30 −20 −10 10 20 30 gdB =−12dB gdB =−6dB gdB =0dB gdB =6dB gdB =12dB Bass EQ: Amplitude Response freqeuncy(log10) Gain(dB) 10

2

10

3

10

4

−30 −20 −10 10 20 30 gdB =−12dB gdB =−6dB gdB =0dB gdB =6dB gdB =12dB Treble EQ: Amplitude Response freqeuncy(log10) Gain(dB)

Bi-quad Filters

§ Low-pass filter

– fc : cut-off frequency, Q: resonance

36

H(z) = (1−cosΘ 2 ) 1+ 2z−1 +1z−2 (1+α)− 2cosΘz−1 +(1−α)z−2 α = sinΘ 2Q

10

2

10

3

10

4

−30 −20 −10 10 20 30 f=400 f=1000 f=3000 f=8000 Lowpass Filters freqeuncy(log10) Gain(dB) 10

2

10

3

10

4

−30 −20 −10 10 20 30 Q =0.5 Q =1 Q =2 Q =4 Lowpass Filters freqeuncy(log10) Gain(dB)

Θ = 2π fc / fs

Bi-quad Filters

§ High-pass filter

37

H(z) = (1+cosΘ 2 ) 1− 2z−1 +1z−2 (1+α)− 2cosΘz−1 +(1−α)z−2 α = sinΘ 2Q Θ = 2π fc / fs

10

2

10

3

10

4

−30 −20 −10 10 20 30 f=400 f=1000 f=3000 f=8000 Highpass Filters freqeuncy(log10) Gain(dB) 10

2

10

3

10

4

−30 −20 −10 10 20 30 Q =0.5 Q =1 Q =2 Q =4 Highpass Filters freqeuncy(log10) Gain(dB)

Bi-quad Filters

§ Band-pass filter

38

H(z) = (sinΘ 2 ) 1− z−2 (1+α)− 2cosΘz−1 +(1−α)z−2

10

2

10

3

10

4

−30 −20 −10 10 20 30 f=400 f=1000 f=3000 f=8000 Bandpass Filters freqeuncy(log10) Gain(dB) 10

2

10

3

10

4

−30 −20 −10 10 20 30 Q =0.5 Q =1 Q =2 Q =4 Bandpass Filters freqeuncy(log10) Gain(dB)

α = sinΘ 2Q Θ = 2π fc / fs

Bi-quad Filters

§ Notch filter

39

H(z) = 1− 2cosΘz−1 + z−2 (1+α)− 2cosΘz−1 +(1−α)z−2 α = sinΘ 2Q Θ = 2π fc / fs

10

2

10

3

10

4

−30 −20 −10 10 20 30 f=400 f=1000 f=3000 f=8000 Notch Filters freqeuncy(log10) Gain(dB) 10

2

10

3

10

4

−30 −20 −10 10 20 30 Q =0.5 Q =1 Q =2 Q =4 Notch Filters freqeuncy(log10) Gain(dB)

10

2

10

3

10

4

−30 −20 −10 10 20 30 AdB=−12 AdB=−6 AdB=0 AdB=6 AdB=12 EQ freqeuncy(log10) Gain(dB) 10

2

10

3

10

4

−30 −20 −10 10 20 30 AdB=−12 AdB=−6 AdB=0 AdB=6 AdB=12 EQ freqeuncy(log10) Gain(dB)

Bi-quad Filters

§ Equalizer

40

H(z) = (1+α ⋅ A)− 2cosΘz−1 +(1+α ⋅ A)z−2 (1+α / A)− 2cosΘz−1 +(1−α / A)z−2 α = sinΘ 2Q Θ = 2π fc / fs

Q=1 Q=4

References

§ Cookbook formulae for audio EQ biquad filter coefficient,

• R. Bristow-Johnson

– http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt

41