**Hilbert transform system:**

The delay of (filter tap length)/2 is given to I-channel of IQ-signals. It and the Hilbert transform output of Q-signal are subtracted or added. The direction separated signals are calculates by formulas (9) and (10). Here a convolution is indicated . The tap number is set to 128 in the estimation of the calculation load shown in Table 3.

$$F1(n) = \operatorname{Re}\left(X(n - n \text{trap } / \, \mathcal{D})\right) + \operatorname{Im}\left(X(n) \otimes h1(n \text{trap})\right) \tag{9}$$

$$R1(n) = -\operatorname{Re}\left(X(n - n\text{trap }/\ \mathcal{D})\right) + \operatorname{Im}\left(X(n) \otimes h1(n\text{trap})\right) \tag{10}$$

The coefficient *h1* of Hilbert transform is given by a formula (11).

Complex Digital Filter Designs for Audio Processing in Doppler Ultrasound System 219

Z

IQ signal is modulated and frequency is shifted *+fs/4* and *–fs/4*. The positive-component (0 to *+fs/2*) and negative-component (*-fs/2* to *0*) are extracted by LPF. The *+fs/4* shift and the *–fs/4* shift are returned by demodulation. The direction separation outputs are calculated by formulas (20) and (21). The example of Table 3 is referred to the prior art. The 128-tap FIR

> 5( ) ( ) exp ( ) exp 2 2 *F n Xn j n CLPF ntap j n* § · § · § · S

> 5( ) ( ) exp ( ) exp 2 2 *R n Xn j n CLPF ntap j n* § · § · § · S

The direction separation outputs are calculated by formulas (22) and (23).

There are two sets of phase-shifter with the transfer characteristic that makes relative phase difference of IQ-signal 90 degree. The addition-and-subtraction of these outputs is used.

*F z X z Phase z X z Phase z* 6( ) Re ( ) 1( ) Im ( ) 2( )   (22)

The two sets of phase-shifter are the cascade connection of second-order all-pass filter arrays. They are denoted by formulas (24) and (25) as a *Phase1 (z)* and a *Phase2 (z).* In the estimation of Table 3, the cascade connections of four steps of all-pass filters are used. Moreover, in order to improve the performance near the Nyquist frequency, an interpolator and a decimator are added before and after phase-shifter. Table 3 is calculated in *N*= 4, and

 (19)

 S

 S

§ · ¨ ¸ ¨ ¸ ¨ ¸ ¨ ¸ © ¹ © ¹ © ¹ © ¹ (20)

§ · ¨ ¸ ¨ ¸ ¨ ¸ ¨ ¸ © ¹ © ¹ © ¹ © ¹ (21)

*R z X z Phase z X z Phase z* 6( ) Re ( ) 1( ) Im ( ) 2( )   (23)

1

1

*<sup>n</sup> <sup>k</sup> k k*

 

*a z*

*b z*

*<sup>n</sup> <sup>k</sup> k k*

 

1 1( ) <sup>1</sup>

1 2( ) <sup>1</sup>

*z a Phase z*

*z b Phase z*

Above six kinds of signal-processing algorithms are confirmed by the simulation. The chirpwaveform that frequency and a direction are changed is used as an input. The result of a simulation is shown in Fig. 7. Fig. 7(a) is an input signal and the sign of frequency has inverted near 200ck (equivalent to the time shown in the Fig. 7 broken line). A solid line is Isignal and a dotted line is Q-signal. Figures 7(b) to (g) are output waveforms of each signalprocessing system. A solid line is a positive-output (forward) of the Doppler audio, and a dotted line is a negative-output (reverse) of the Doppler audio. Amplitude of positiveoutput becomes large on the right-hand side of a broken line, and it becomes small on the

1

1

(24)

(25)

*R n IFFT WR FFT X n h n* 4( ) Re ( ) ( ) 4( )

**Modulation/demodulation system:** 

**Phase-shift system:** 

low-pass filter, which has 63/128 cut-off, is used.

the FIR filter of *2N* tap is used as an interpolator.

$$\begin{aligned} h1(n) &= \frac{2}{\pi} \cdot \frac{\sin^2(\pi \cdot n / 2)}{n} \quad &\text{( $n \neq 0$ )}\\ &= 0 \quad &\text{( $n = 0$ )} \end{aligned} \tag{11}$$

#### **Complex FIR system:**

There is a report of the Doppler audio separation processing using a complex FIR filter. However, since there is no description about a filter coefficient, we designed in a frequency domain and transformed into FIR coefficient in time domain using inverse Fourier transform. The output of complex FIR system is denoted by formulas (12) and (13). In the estimation of Table 3, the 128-tap coefficient sequence with the pass band of r *fs* /128 to r 63 /128 *fs* is used.

$$F2(n) = X(n) \otimes HF2(n \text{top}) \tag{12}$$

$$R\mathfrak{D}(n) = X(n) \otimes HR\mathfrak{D}(n \text{top}) \tag{13}$$

#### **Complex IIR system:**

Based on the shift theory of Fourier transform, frequency shift is applied to *z* operators. A real-LPF transfer function is changed into the positive-BPF and the negative-BPF. The complex IIR transfer functions become a formulas (14) and (15).

$$F\mathfrak{Z}(z) = HF\mathfrak{Z}(z) \cdot X(z) \tag{14}$$

$$R\mathfrak{B}(z) = HR\mathfrak{B}(z) \cdot X(z) \tag{15}$$

When the transfer function of real LPF is set to *RLPF (z),* transfer functions of *HF3 (z')* and *FR3 (z'')* are calculated by transformed operators. In the estimation of Table 3, the filter with the 8th order Butterworth type is used.

$$HF3(z) = RLPFS(z') \qquad \text{where} \qquad z' = -j \cdot z \tag{16}$$

$$\text{RHS}(z) = \text{RLPF3}(z'') \quad \text{where} \quad z'' = \mathbf{j} \cdot \mathbf{z} \tag{17}$$

#### **FFT/IFFT system:**

The IQ-signal is separated by the positive-filter and negative-filter in a frequency domain. Next, the separated spectra are returned to waveforms in time domain by inverse-FFT. There is a report of this system aiming at the Doppler noise rejection. For the continuous output after inverse-FFT, a shift addition of the time waveform is carried out in time domain. The outputs of this system can be denoted by formulas (18) and (19). In estimation of Table 3, FFT/IFFT point number is set to 128, and used the frequency filter of r *fs* /128 to r 63 /128 *fs* for separation. Moreover, Hamming window (*h4*) is applied, and 32 timeseries are shift-added.

$$F4(n) = \operatorname{Re}\left\{IFFT\left(VNF(o) \cdot FFT\left(X(n)\right)\right)\right\} \cdot h4(n) \tag{18}$$

$$R4(n) = \operatorname{Re}\left(IFFT\left(VRR(a) \cdot FFT\left(X(n)\right)\right)\right) \cdot h4(n) \tag{19}$$

#### **Modulation/demodulation system:**

IQ signal is modulated and frequency is shifted *+fs/4* and *–fs/4*. The positive-component (0 to *+fs/2*) and negative-component (*-fs/2* to *0*) are extracted by LPF. The *+fs/4* shift and the *–fs/4* shift are returned by demodulation. The direction separation outputs are calculated by formulas (20) and (21). The example of Table 3 is referred to the prior art. The 128-tap FIR low-pass filter, which has 63/128 cut-off, is used.

$$F5(n) = \left( \left( X(n) \cdot \exp\left(-\frac{\pi}{2} \cdot j \cdot n\right) \right) \otimes \text{CLPF}(n \text{tap}) \right) \cdot \exp\left(+\frac{\pi}{2} \cdot j \cdot n\right) \tag{20}$$

$$R\mathbf{5}(n) = \left( \left( X(n) \cdot \exp\left( + \frac{\pi}{2} \cdot j \cdot n \right) \right) \otimes \text{CLPF}(n \text{atap}) \right) \cdot \exp\left( -\frac{\pi}{2} \cdot j \cdot n \right) \tag{21}$$

#### **Phase-shift system:**

218 Applications of Digital Signal Processing

*<sup>n</sup> h n <sup>n</sup> n*

There is a report of the Doppler audio separation processing using a complex FIR filter. However, since there is no description about a filter coefficient, we designed in a frequency domain and transformed into FIR coefficient in time domain using inverse Fourier transform. The output of complex FIR system is denoted by formulas (12) and (13). In the estimation of Table 3, the 128-tap coefficient sequence with the pass band of r *fs* /128 to

Based on the shift theory of Fourier transform, frequency shift is applied to *z* operators. A real-LPF transfer function is changed into the positive-BPF and the negative-BPF. The

When the transfer function of real LPF is set to *RLPF (z),* transfer functions of *HF3 (z')* and *FR3 (z'')* are calculated by transformed operators. In the estimation of Table 3, the filter with

*HF z RLPF z* 3( ) 3( ') where *z jz* ' (16)

The IQ-signal is separated by the positive-filter and negative-filter in a frequency domain. Next, the separated spectra are returned to waveforms in time domain by inverse-FFT. There is a report of this system aiming at the Doppler noise rejection. For the continuous output after inverse-FFT, a shift addition of the time waveform is carried out in time domain. The outputs of this system can be denoted by formulas (18) and (19). In estimation of Table 3, FFT/IFFT point number is set to 128, and used the frequency filter of r *fs* /128 to r 63 /128 *fs* for separation. Moreover, Hamming window (*h4*) is applied, and 32 time-

> *F n IFFT WF FFT X n h n* 4( ) Re ( ) ( ) 4( )  Z

complex IIR transfer functions become a formulas (14) and (15).

S

**Complex FIR system:** 

r 63 /128 *fs* is used.

**Complex IIR system:** 

**FFT/IFFT system:** 

series are shift-added.

the 8th order Butterworth type is used.

<sup>2</sup> 2 sin ( / 2) 1( ) ( 0)

z

S

0 ( 0)

*n*

*F n X n HF ntap* 2( ) ( ) 2( ) (12)

*R n X n HR ntap* 2( ) ( ) 2( ) (13)

*F z HF z X z* 3( ) 3( ) ( ) (14)

*R z HR z X z* 3( ) 3( ) ( ) (15)

  (18)

*HR z RLPF z* 3( ) 3( '') where *z jz* '' (17)

(11)

There are two sets of phase-shifter with the transfer characteristic that makes relative phase difference of IQ-signal 90 degree. The addition-and-subtraction of these outputs is used. The direction separation outputs are calculated by formulas (22) and (23).

$$F\mathfrak{G}(z) = \operatorname{Re}\left(X(z)\right) \cdot \operatorname{Phase1}(z) + \operatorname{Im}\left(X(z)\right) \cdot \operatorname{Phase2}(z) \tag{22}$$

$$R\Theta(z) = -\operatorname{Re}\left(X(z)\right) \cdot \operatorname{Phase1}(z) + \operatorname{Im}\left(X(z)\right) \cdot \operatorname{Phase2}(z) \tag{23}$$

The two sets of phase-shifter are the cascade connection of second-order all-pass filter arrays. They are denoted by formulas (24) and (25) as a *Phase1 (z)* and a *Phase2 (z).* In the estimation of Table 3, the cascade connections of four steps of all-pass filters are used. Moreover, in order to improve the performance near the Nyquist frequency, an interpolator and a decimator are added before and after phase-shifter. Table 3 is calculated in *N*= 4, and the FIR filter of *2N* tap is used as an interpolator.

$$Phase1(z) = \prod\_{k=1}^{n} \frac{z^{-1} - a\_k}{1 - a\_k \cdot z^{-1}} \tag{24}$$

$$Phase2(z) = \prod\_{k=1}^{n} \frac{z^{-1} - b\_k}{1 - b\_k \cdot z^{-1}} \tag{25}$$

Above six kinds of signal-processing algorithms are confirmed by the simulation. The chirpwaveform that frequency and a direction are changed is used as an input. The result of a simulation is shown in Fig. 7. Fig. 7(a) is an input signal and the sign of frequency has inverted near 200ck (equivalent to the time shown in the Fig. 7 broken line). A solid line is Isignal and a dotted line is Q-signal. Figures 7(b) to (g) are output waveforms of each signalprocessing system. A solid line is a positive-output (forward) of the Doppler audio, and a dotted line is a negative-output (reverse) of the Doppler audio. Amplitude of positiveoutput becomes large on the right-hand side of a broken line, and it becomes small on the

Complex Digital Filter Designs for Audio Processing in Doppler Ultrasound System 221

The comparison of time-delay is shown in Table 2. Frequency resolution is adjusted by parameter of each system in accordance with the target performance of Table 1. Since the signal-processing inputs are sampled by *fs*, time-delay will become large if *fs* becomes low. Table 2 is calculated by *fs*=4kHz condition. Incidentally by *fs*=1kHz, time-delay increases 4 times. The time-delay caused by operation is assumed zero, and estimated only the delay caused by sampling simply. Moreover, the influence of the transient response of the complex IIR system and the phase-shift system is not taking into

(a) input signal

forward

reverse

Fig. 8. Comparison of response between complex FIR method and complex IIR method

method estimation time-delay (ms) Hilbert transform *tap/fs* 32 (*tap*=128) Complex FIR *tap/fs* 32 (*tap*=128) Complex IIR *order/fs* (\*1) 2 (*order*=8) FFT/IFFT *1.5\*N/fs* (\*2) 48 (*N*=128) Moduration/Demoduration *tap/fs* 32 (*tap*=128) Phase shift *max(2,order/N)/fs* (\*1) 1 (*order*=4, *N*=1) Estimated at *fs*=4kHz, (\*1) not including transient response, (\*2) IFFT shift addition pitch is *N/4*

I-channel Q-channel

(c) output signal of the complex IIR system

As calculation load depends on the hardware architecture, the multiplication and addition times per 1 second (floating point single precision) is used for this estimation. Moreover, the complex-multiplication is considered as 4 times, and complex-addition is considered as twice. The overhead of the processing which requires a lot of memory buffers is assumed to be 20%. The other overhead is assumed to be 10%. The calculation elements, estimation formula and calculation load for each signal-processing systems are shown in Table 3. Incidentally, at *fs*=52 kHz (maximum PRF in actual system), calculation load increases 13 times. The result of Table 2 and Table 3 shows that the complex IIR system and the phase-

(b) output signal of the complex FIR system

reverse

Time (ck)

Time (ck)

Time (ck)

consideration here.

Amplitude

Amplitude

Amplitude

Table 2. Comparison of time-delay

forward

left-hand side of the broken line. Amplitude of negative-output becomes small on the righthand side of a broken line, and it becomes large on the left-hand side of the broken line. This result shows that each system works correctly. Moreover, it shows that the waveform and response time at the turning point of sign (near the DC) have a difference among the systems. As these causes, performance differences, such as the response characteristic and delay time, can be considered.

Fig. 7. Chirp wave responses

### **3.3 Comparison of time-delay and calculation load**

The response is important for blood vessel detection, and the time-delay estimates it. The simulation result of the time-delay is shown in Fig. 8. They are response waveforms of the sinusoidal-waveform that changes discontinuously. Solid line and dotted line are I-signal and Q-signal in Fig. 8(a). Amplitude and frequency are changing near the 50ck. The solid lines of Fig. 8(b) and Fig. 8(c) are positive-output waveforms, and dotted lines are negativeoutput waveforms. The output waveform of the complex FIR system of Fig. 8 (b) changed from a turning point of the input shown with the dashed line after 64ck (time shown with the chain line among Fig. 8(b)), and is stable gradually. The output waveform of the complex IIR system of Fig. 8(c) is stable from the turning point after 8ck (time shown with the chain line among Fig. 8(c)).

220 Applications of Digital Signal Processing

left-hand side of the broken line. Amplitude of negative-output becomes small on the righthand side of a broken line, and it becomes large on the left-hand side of the broken line. This result shows that each system works correctly. Moreover, it shows that the waveform and response time at the turning point of sign (near the DC) have a difference among the systems. As these causes, performance differences, such as the response characteristic and

I-channel

Q-channel

(b) output signal of the Hilbert transform system

reverse

forward

reverse

forward

(c) output signal of the complex FIR system

forward

reverse

(d) output signal of the complex IIR system

(e) output signal of the FFT/IFFT system

forward

reverse

(g) output signal of the phase-shift system

The response is important for blood vessel detection, and the time-delay estimates it. The simulation result of the time-delay is shown in Fig. 8. They are response waveforms of the sinusoidal-waveform that changes discontinuously. Solid line and dotted line are I-signal and Q-signal in Fig. 8(a). Amplitude and frequency are changing near the 50ck. The solid lines of Fig. 8(b) and Fig. 8(c) are positive-output waveforms, and dotted lines are negativeoutput waveforms. The output waveform of the complex FIR system of Fig. 8 (b) changed from a turning point of the input shown with the dashed line after 64ck (time shown with the chain line among Fig. 8(b)), and is stable gradually. The output waveform of the complex IIR system of Fig. 8(c) is stable from the turning point after 8ck (time shown with

(f) output signal of the modulation/demodulation system

reverse

forward

forward

reverse

Time (ck)

(a) input signal

delay time, can be considered.

Amplitude

Amplitude

Amplitude

Amplitude

Amplitude

Amplitude

Amplitude

**3.3 Comparison of time-delay and calculation load** 

Fig. 7. Chirp wave responses

the chain line among Fig. 8(c)).

The comparison of time-delay is shown in Table 2. Frequency resolution is adjusted by parameter of each system in accordance with the target performance of Table 1. Since the signal-processing inputs are sampled by *fs*, time-delay will become large if *fs* becomes low. Table 2 is calculated by *fs*=4kHz condition. Incidentally by *fs*=1kHz, time-delay increases 4 times. The time-delay caused by operation is assumed zero, and estimated only the delay caused by sampling simply. Moreover, the influence of the transient response of the complex IIR system and the phase-shift system is not taking into consideration here.

Fig. 8. Comparison of response between complex FIR method and complex IIR method


Estimated at *fs*=4kHz, (\*1) not including transient response, (\*2) IFFT shift addition pitch is *N/4*

### Table 2. Comparison of time-delay

As calculation load depends on the hardware architecture, the multiplication and addition times per 1 second (floating point single precision) is used for this estimation. Moreover, the complex-multiplication is considered as 4 times, and complex-addition is considered as twice. The overhead of the processing which requires a lot of memory buffers is assumed to be 20%. The other overhead is assumed to be 10%. The calculation elements, estimation formula and calculation load for each signal-processing systems are shown in Table 3. Incidentally, at *fs*=52 kHz (maximum PRF in actual system), calculation load increases 13 times. The result of Table 2 and Table 3 shows that the complex IIR system and the phase-

Complex Digital Filter Designs for Audio Processing in Doppler Ultrasound System 223

separation performance. The separation performance is deteriorated especially near the

*tap* =128 *tap* =32 *tap* =8

Power (dB)

forward forward

Power (dB)

forward forward

Power (dB)

forward forward

Fig. 9. Example of frequency response: the Hilbert transform system

(a) the Hilbert transform system

Frequency (*fs*) (c) the complex IIR system

(e) the modulation/demodulation

Frequency (*fs*)

reverse reverse

reverse reverse

Frequency (*fs*)

We made the target performances of the direction separating process of digital Doppler audio, and evaluated six kinds of digital-signal-processing ideas that were pre-existing or were newly devised. The performances of each processing were evaluated by comparing

1. The complex IIR system and the phase-shift system are filling the target performance of

2. The target performance of direction separation is filled except for the phase-shift

Fig. 10. Frequency characterization and direction separation performance

many responses such as chirp or step and so on. The results are following.

reverse reverse

Frequency (*fs*)

(b) the complex FIR system

(d) the FFT/IFFT system

(f) the phase-shift system

Frequency (*fs*)

Frequency (*fs*)

Frequency (*fs*)

Nyquist frequency.

Power (dB)

Power (dB)

Power (dB)

**3.5 Conclusion** 

response time.

system.

Power (dB)


shift system are filling the target performance of time-delay. It turns out that calculation load is light in order of the phase-shift system, the complex IIR system, and the Hilbert transform system.

R-add: real addition, R-mul: real multiplication, Ovh: over head, C-add: complex addition, C-mul: complex multiplication, Calculation load is estimated at *fs*=4kHz

Table 3. Comparison of calculation load
