Identification and localization of the sources of cyclostationary sound fields

Applied Acoustics 87 (2015) 64–71

Contents lists available at ScienceDirect

Applied Acoustics journal homepage: www.elsevier.com/locate/apacoust

Identification and localization of the sources of cyclostationary sound fields Zhi Min Chen a,b,⇑, Hai Chao Zhu a,b, Min Peng c a

Institute of Noise & Vibration, Naval University of Engineering, Wuhan 430033, PR China National Key Laboratory on Ship Vibration & Noise, Wuhan 430033, PR China c Engineering Institute, Wuhan Yangtze Business University, Wuhan 430065, PR China b

a r t i c l e

i n f o

Article history: Received 30 April 2013 Received in revised form 31 March 2014 Accepted 19 June 2014

Keywords: Cyclostationary sound field Planar near-field acoustic holography The cyclic spectral density Identification and localization of the sources

a b s t r a c t A cyclostationary sound field is approximately regarded as a stationary sound field when analyzed by the traditional planar near-field acoustic holography (PNAH), which ignores the periodical time-variant property and makes it hard to reflect accurately the radiation characteristics based on the hologram results. In this article, a cyclostationary planar near-field acoustic holography (CPNAH) technique is proposed in which the cyclic spectral density (CSD) instead of the complex sound pressure is adopted as reconstructing physical quantity and the physical properties of CSD are utilized. And to avoid plenty of calculation and improve accuracy of extracting the cyclic property, pre-processing of holographic data, the method known as the gathering slices of CSD is also proposed. The experiment results demonstrate that the cyclic properties of cyclostationary sound fields can be extracted and the sources of cyclostationary sound fields can be exactly identified and localized by means of CPNAH. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction The cyclostationary sound field is a special kind of non-stationary sound field which exists widely in engineering practice. For example, the sound field radiated from the rotating machinery or the reciprocating machinery has a significant periodical timevariant characteristic due to the symmetrical or approximate symmetrical structure and the periodic operating mode of the machine [1–4]. Such a sound field is usually treated as a stationary sound field and its properties in frequency domain are analyzed by the Fourier analysis tool. This processing method is simple by which some engineering problems have been solved, however, the loss of some information about the frequency change with the passage of time cannot be avoided. The cyclic statistics theory is an appropriate analysis tool for cyclostationary sound fields [5,6]. The non-stationary character i.e. periodical time-variant properties extracted by means of analyzing the cyclic statistics of the cyclostationary sound field can fully reflect the operating state and physical properties of the rotating machinery or the reciprocating machinery. Nowadays the second-order cyclic statistics have been used in engineering ⇑ Corresponding author at: Institute of Noise & vibration, Naval University of Engineering, Wuhan 430033, PR China. Tel.: +86 27 83443247; fax: +86 27 83443981. E-mail address: [email protected] (Z.M. Chen). http://dx.doi.org/10.1016/j.apacoust.2014.06.013 0003-682X/Ó 2014 Elsevier Ltd. All rights reserved.

practice more widely than other cyclic statistics, and some tangible results have been obtained. In near-field acoustic holography (NAH), an acoustic quantity in near field, such as complex sound pressure, is used to reconstruct the sound field on the surface of the radiator for the sake of identifying and localizing the sound source [7,8]. However the sound field is usually treated as the stationary field to this day when it is generated by the rotating machinery or the reciprocating machinery, this is restricted to some extent. A near-field acoustic holography technique for cyclostationary sound fields was presented in Ref. [9]. In the research, a pair of cyclic cross-spectral density functions of the sound pressure at holographic measuring points and at reference point were adopted as the reconstruction physical quantity of NAH and two times space reconstructions were carried out to get the CSD distribution on the reconstruction plane. The traditional time domain smooth periodogram method was used for the pre-processing of measured signals. Obviously, it is difficult to adopt the complicated technique. The CSD, which belongs to the second-order cyclic statistics, is adopted as the reconstruction quantity of NAH in this article. The advantage is that the cyclic components of cyclostationary sound field may be extracted accurately, noise suppressed greatly, and phase information of the signal reserved [10,11]. The disadvantage lies in its calculation complexity and a long length of data accumulation in order to obtain higher estimation accuracy. To

65

Z.M. Chen et al. / Applied Acoustics 87 (2015) 64–71

2.2. PNAH of cyclostationary sound field In PNAH technique, the data on hologram plane is sampled within a two-dimension aperture by scanning and the holography reconstruction is actually a problem of data processing based on two-dimensional matrix. When the complex sound pressure is replaced by the CSD as the reconstruction physical quantity, according to the physical character of CSD, the delay effect of the measurement data between two close scanning steps must be eliminated to ensure spatial phase relationship of the two-dimensional holography data matrix. The power spectral density (PSD) is insensitive to time delay, whereas the CSD is sensitive due to the presence of cyclic frequency domain. Suppose a signal is v(t) = x(t  sn), where sn is the time delay. The CSD of v(t) is

Sav ðf Þ ¼ Sax ðf Þej2pasn Fig. 1. Amplitude of CSD.

solve the calculation problem, a gathering slices method of CSD is proposed by referring to time aliasing methods on time series [12]. With the help of analyzing zero delay slice of the cyclic autocorrelation function of the signal, the peak frequencies of interest on the cyclic frequency axis are selected to calculate CSD slices, and the interference between overlapping terms is eliminated by analyzing gathering CSD slices, consequently feature frequencies of sound signal are extracted exactly and the planar NAH (PNAH) reconstruction based on two-dimensional Fourier transform is carried out at feature frequencies. Generally speaking, the advantages of both PNAH and CSD are employed in this method, and feature frequencies of the sound field can be extracted, while the sound energy at feature frequencies can be visualized. 2. Theory 2.1. PNAH of stationary sound field Suppose that the complex sound pressure on hologram plane is p(xh, yh, zh) and Green’s function under Dirichlet boundary conditions is GD(x, y, z), the complex sound pressure on any plane parallel with hologram plane can be expressed as

pðx; y; zÞ ¼

ZZ s

ðpðxh ; yh ; zh ÞGD ðx  xh ; y  yh ; z  zh ÞÞdS

ð1Þ

from which it can be seen that a delay factor e in cyclic frequency domain (a domain) exists between the initial signal and the delayed signal. The delay factor between two close scanning steps can be obtained using a fixed reference microphone in scanning process. Suppose the pressure measured by the reference microphone at the nth scanning step is ^rn ðtÞ ¼ rðt  sn Þ, where n = 1:N, s1 = 0, N is total step number. The CSD of ^rn ðtÞ is S^ran ðf Þ, then the delay factor can be obtained by the following expression,

ej2pasn ¼ S^ran ðf Þ=S^ra1 ðf Þ

Sapn ðf Þ ¼ Spa^n ðf ÞS^ra1 ðf Þ=S^arn ðf Þ

e D ðkx ; ky ; z  zh Þ ~ðkx ; ky ; zÞ ¼ F 1 ½p ~ðkx ; ky ; zh Þ G pðx; y; zÞ ¼ F 1 ½p

0.6

0.7

Once the complex sound pressure on a near field plane is obtained by scanning, the complex sound pressure on any plane parallel to the near field plane can be obtained by Eq. (2).

Rxα (0) /um

2

ð2Þ

ð3Þ

ð6Þ

When time delay between two close scanning steps is eliminated by Eq. (6), measuring data on the hologram plane is equivalent to synchronous acquisition, and phases of CSD at the cyclic frequencies are reserved, then the spatial phases of data matrices on the two-dimensional hologram plane exist. So if the CSD rather than the complex sound pressure is adopted as the reconstruction variable, the reconstruction can be carried out directly using Eq. (2). This technique is named cyclostationary planar near-field acoustic holography (CPNAH). There exist some substantial differences between the PNAH technique and the CPNAH technique in

0.8

e D ðkx ; ky ; z  zh Þ express ~ðkx ; ky ; zÞ, p ~ðkx ; ky ; zh Þ and G where p two-dimensional Fourier transform of p(x, y, z), p(xh, yh, zh) and GD(x, y, z  zh) respectively. Eq. (2) is the space transformation formula of the sound field. e D ðkx ; ky ; z  zh Þ may be derived by function integral table or the G wave equation, the result is

ð5Þ

Let the CSD of signals measured by the microphone array at the nth scanning step be Spa^n ðf Þ. After eliminating the effect of time delay, the CSD can be written as

The two-dimensional Fourier transform is carried out on both sides of Eq. (1) and the convolution theorem is used, the following expressed can be derived,

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 8 < ejðzzh Þ k2 ðk2x þk2y Þ ; k2  k2 þ k2 x y e D ðkx ; ky ; z  zh Þ ¼ G pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : ðzzh Þ ðk2x þk2y Þk2 2 2 2 e ; k < kx þ ky

ð4Þ j2pasn

0.5 0.4 0.3 0.2 0.1 0 -1000

-500

0

500

α/Hz Fig. 2. The cyclic autocorrelation function (l = 0).

1000

66

Z.M. Chen et al. / Applied Acoustics 87 (2015) 64–71

xðaÞ ðnÞ ¼

M X xðn þ mNÞ n ¼ 0; 1; 2; . . . ; N  1

ð7Þ

m¼0

The time aliasing process in the time domain can effectively improve the computational efficiency and the accuracy of spectrum estimation. Therefore, the time aliasing process is conducted on the time series of time-variant autocorrelation function. The non-symmetric time-variant autocorrelation function of x(n) is

Rx ½n; l ¼ x½nx½n þ l

ð8Þ

Time aliasing process is carried out on each time series with fixed time delay l in Eq. (8), then the time aliasing signal standardized by 1/M of the time-variant autocorrelation function is

RðaÞ x ½n; l ¼

1 1 X X 1M 1M Rx ½n þ mN; l ¼ x½n þ mNx½n þ mN þ l M m¼0 M m¼0

ð9Þ

Fig. 3. The representation of the gathering slices of CSD.

After the Fourier transform being carried out on the time variable of RðaÞ x ½n; l, the cyclic autocorrelation function (CAF) is which the reference microphone is used to eliminate the time delay and the reconstruction result accessible is the CSD of sound signals. 3. Pre-processing of data measured on the holography plane The key to CPNAH technique is the pre-processing of data measured on the holography plane so that spectral frequencies indicating modulation characteristics of the sound field may be extracted accurately and effectively by CSD. Most of existing fast algorithms for CSD, such as the fast Fourier transform accumulation method (FAM) [13], may be classed as periodogram methods similar to the calculation of power spectral density. However, the computation of such algorithms is still complex, and cyclic frequencies and spectral frequencies of interest cannot be selected directly from three-dimensional map of CSD, especially there exist nonzero components at the non-cyclic frequencies due to the limited length of the measured data. Therefore it is difficult to extract feature frequencies. To solve this problem, a gathering slices method of CSD is presented. Suppose a discrete cyclostationary signal is x(n), its time aliasing signal after time aliasing process in the time domain is x(a)(n). The time aliasing process is a repeating process with period N samples that has been created by summing together sample points from M successive, N point blocks of the original sampled signal, and may be described mathematically as

Rax ðlÞ ¼

N1 X j2panT s jpalT s RðaÞ e x ½n; le

ð10Þ

n¼0

where a = n/NTs is the cyclic frequency, Ts is the sampling interval, and ejpalT s is the compensation for shifting of non-symmetric frequency. The Fourier transform is carried out on the time delay variable l in Eq. (10), the CSD is

Sax ðf Þ ¼

w1 X Rax ðlÞej2pflT s

ð11Þ

l¼0

in which w is the width of the time delay window. It can be found from Eqs. (10) and (11) that the CAF and the CSD of the signal has the same cyclic frequency set. Rax ðlÞ takes the maximum value in a domain when l = 0. Thus, the cyclic frequency set can be obtained by the zero delay slices of CAF of the signal, and cyclic frequencies of interest may be chosen to calculate CSD slices. The analyzing process using the gathering slices method of CSD is summarized as follows: (1) To calculate zero delay slices of cyclic autocorrelation function (l = 0) using Eq. (10), and select peak cyclic frequencies of interest to form a cyclic frequency set. (2) To substitute a cyclic frequency in the cyclic frequency set into Eq. (10) to obtain the slice of cyclic autocorrelation function at this cyclic frequency. (3) To perform Fourier transform in the time delay domain on each slice of CAF to obtain the gathering slices of CSD. To prove the validity of the gathering slices method of CSD, let us conduct a numerical simulation on a typical cyclostationary signal. Assume an amplitude modulated signal is

xðtÞ ¼ ð1 þ cosð2pf2 tÞÞ cosð2pf1 tÞ

Fig. 4. The scene of experiment.

ð12Þ

where the carrier frequency f1 = 300 Hz, amplitude modulating frequency f2 = 32 Hz. Obviously this is a typical cyclostationary signal. The sampling frequency is 4096 Hz and the sampling time is 10 s. The numerical simulation results are shown in Figs. 1–3. The CSD is calculated using the periodogram method of timedomain smoothing algorithm firstly. Its contour map is shown in Fig. 1. The support domain of CSD in the dual-frequency (a, f) plane is a diamond-shaped area, and the CSD at non-cyclic frequencies is zero. The zero delay slices of CAF of the amplitude modulated signal are analyzed, the result is shown in Fig. 2. Sax ðf Þ is calculated

67

Z.M. Chen et al. / Applied Acoustics 87 (2015) 64–71

(b) -100

(a) -100

-150

-200

PSD/Pa2.Hz-1

PSD/Pa2.Hz-1

-150

-250 -300

-250

-300

-350 -400

-200

0

100 200 300 400 500 600 700 800 900 1000

f/Hz

-350

0

100 200 300 400

500 600 700

800 900 1000

f/Hz

Fig. 5. Power spectral density of reference signal: (a) the first case and (b) the second case.

respectively at each peak cyclic frequency a = 0, 32 Hz, 64 Hz, 536 Hz, 568 Hz, 600 Hz, 632 Hz, 664 Hz, the gathering slices are shown in Fig. 3. The information implied in Fig. 3 is identical to that in Fig. 1, but plenty of calculation of spectral correlation at a large number of non-cyclic frequencies is avoided, and the cyclic frequencies and spectral frequencies can be selected directly by means of the gathering slices method of CSD. The method is superior to periodogram methods in reducing calculating complexity, feature extracting accuracy and intuitive graph. When a = 0, the CSD is degraded into the traditional power spectral density and modulating frequency exhibits the sideband of carrier frequency. 4. Experimental investigations 4.1. Experiment on a motor The experiment was conducted on an Y280M-6 type threephase asynchronous motor fixed on an isolated platform. The rated speed of the motor is 980 r/min. A self-developed measuring system for NAH was used to sample data by scanning method. The measuring grid comprised 19 (horizontal direction)  12 (vertical direction) points with lattice spacing of 8 cm in both directions. The holography plane was 33 cm away from the surface of the

motor. To ensure the accuracy of reconstruction, the holographic plane distance from the source to smaller than a third of wavelength of the analyzing frequency, space sampling point spacing is smaller than the 1/2 wavelength of the analyzing frequency, the holographic measurement parameter settings meet the requirements. The reference microphone was placed in front of the motor. The sampling frequency was 2048 Hz and the sampling time was 10s. The background noise is 50 dB. The experimental scene is shown in Fig. 4. The experiment was conducted in two different cases. In the first case, the motor operated at rated speed with no load. In the second case, the motor operated the same as in the first case, but an iron piece was placed in the way that one end of the piece was fixed to the supporter while the other end was pressed on the surface of the coupling to generate noise by friction (oval area in the experimental scene). The acoustic noise radiated from an electric motor can be categorized mainly into three components: mechanical noise, electromagnetic noise and airborne noise [14]. Mechanical noise is generated mainly by the assembly error between the rotor and bearings. Electromagnetic noise is mainly generated by alternating magnetic forces acting on the motor housing. Airborne noise is mainly generated by the ventilating fan.

Fig. 6. The cyclic autocorrelation function (l = 0): (a) the first case and (b) the second case.

68

Z.M. Chen et al. / Applied Acoustics 87 (2015) 64–71

0.09 0.08 0.07

0.05 0.04 0.03 0.02 0.01 0 -1000 -800 -600 -400 -200

0

200

400

600 800 1000

f/Hz Fig. 7. The CSD slice (a = 16 Hz).

For the sake of comparison, the traditional method for processing stationary sound field was applied firstly. The power spectral density of the sound signal sampled by the reference microphone is shown in Fig. 5 where no obvious change in the second case can be observed by comparing with the first case. It means that the frequency components of the friction noise caused by the iron

-3

0.4 0.3

Y/m

0.4

-0.1

0.1 -0.4

-0.2

0

0.2

0.4

3 2

-0.3

1

-0.4

-0.4 -0.6

4

-0.2

0.2

-0.3

5

0 -0.1

0.3

-0.2

6

0.1

0.5

0

7

0.2

0.6

0.1

8

0.3

0.7

0.2

x 10 9

(b) 0.4

0.8

Y/m

(a)

0.9

-0.6

0.6

-0.4

-0.2

X/m

0

0.2

0.4

0.6

X/m -5

(c)

x 10 14

0.4

x 10 8

(d) 0.4

7

0.3

12

0.3

0.2

10

0.2

6

0.1

5

0

4

-0.1

3

0.1

Y/m

-5

8

0

6

-0.1 -0.2 -0.3

Y/m

S/Pa2.Hz-1

0.06

piece cannot be identified, and of course sound sources cannot be localized accurately. Since the sound field radiated by the motor is time-variant, zero time delay CAF of the reference signal was analyzed, as shown in Fig. 6. It is clear that rotating frequency 16 Hz of the motor appears in the second case. Although frequencies of the friction noise are unknown which are determined by the structure of the iron piece and the interaction between the iron piece and the coupling, the frequency components related to rotating period of the motor would be inherent in the friction noise. Single cyclic spectral density slice at a = 16 Hz was analyzed. The result is shown in Fig. 7 where it can be seen that the peak spectral lines with frequency spacing of 16 Hz appear on the spectral frequency axis. It means that cyclic component of 16 Hz does exist in the slice map. Frequency 188 Hz corresponding to maximum peak was chosen as the analysis frequency to reconstruct the cyclostationary sound field. Let reconstruction frequency be a = 16 Hz and f = 188 Hz, and the reconstruction distance be d = Zc  Zh, where Zc and Zh represent location of the reconstruction plane and the holography plane respectively. The motor center was chosen as the origin and the direction to the holography plane as the positive direction, and two forward and two backward reconstructions at d = ±0.2m, ±0.1m were implemented respectively. The results are displayed in Fig. 8 from which can be seen that the friction noise generated by the iron piece is clearly separated (the left area in the figure). The source in the middle of the motor is generated by electromagnetic

4

-0.2

2

-0.3

2 1

-0.4

-0.4 -0.6

-0.4

-0.2

0

X/m

0.2

0.4

0.6

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

X/m

Fig. 8. The amplitude distribution of CSD (a = 16 Hz, f = 188 Hz): (a) d = 0.2 m, (b) d = 0.1 m, (c) d = 0.1 m and (d) d = 0.2 m.

69

Z.M. Chen et al. / Applied Acoustics 87 (2015) 64–71 -5

0.12

(a) 0.4 0.3

(b) 0.4

0.1

0.2

Y/m

0

0.06

-0.1

0.04

-0.2

0.02

-0.3 -0.4

Y/m

0.08

0.1

x 10 2.2 2

0.3

1.8

0.2

1.6

0.1

1.4 1.2

0

1

-0.1

0.8

-0.2

0.6

-0.3

0.4 0.2

-0.4 -0.6

-0.4

-0.2

0

0.2

0.4

0.6

X/m

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

X/m

Fig. 9. The amplitude distribution of PSD (a = 0, f = 188 Hz): (a) d = 0.2 m and (b) d = 0.2 m.

Fig. 12. The representation of the CSD gathering slice.

Fig. 10. The scene of experiment.

Fig. 11. The cyclic autocorrelation function (l = 0).

excitation acting on the stator and the shell of the motor, and the right source by the ventilating fan. The electromagnetic noise and the airborne noise dominate the radiated noise of the motor at the reconstruction frequency.

Fig. 13. The amplitude distribution of PSD (a = 0, f = 17 Hz).

If the sound field is reconstructed using power spectrum, the periodical components of friction noise will be ignored. Furthermore, because the frequencies of friction noise mix with the rotating frequency, the accurate energy distribution of friction noise

70

Z.M. Chen et al. / Applied Acoustics 87 (2015) 64–71 6

x 10 11

Y/m

0.8

10

0.6

9

0.4

8

0.2

7 6

0

5

-0.2

4

-0.4

3

-0.6

2

-0.8

1 -1

-0.5

0

0.5

1

X/m

Fig. 14. The amplitude distribution of PSD (a = 0, f = 204 Hz).

Fig. 17. The amplitude distribution of CSD (a = 219 Hz, f = 60 Hz).

5

x 10 10

0.8

9

0.6

0.8

8 7

0.2

6

0

5

-0.2

4

-0.4

3

-0.2

-0.6

2

-0.4

1

-0.6

-0.8 -1

-0.5

0

0.5

8

0.6

0.4

7

0.4

6

0.2

Y/m

Y/m

7

x 10 9

5

0

4 3 2 1

-0.8

1

X/m

-1

-0.5

0

0.5

1

X/m

Fig. 15. The amplitude distribution of CSD (a = 17 Hz, f = 25 Hz).

Fig. 18. The amplitude distribution of CSD (a = 407 Hz, f = 237 Hz).

7

x 10 14

x 10 2.2

0.8

0.8 12

0.6 0.4

Y/m

0.2

8

0 6

-0.2 -0.4

4

-0.6

2

-0.8 -1

-0.5

0

0.5

1

X/m Fig. 16. The amplitude distribution of CSD (a = 68 Hz, f = 17 Hz).

Y/m

10

2

0.6

1.8

0.4

1.6

0.2

1.4 1.2

0

1

-0.2

0.8

-0.4

0.6

-0.6

0.4

-0.8

0.2 -1

-0.5

0

0.5

1

X/m Fig. 19. The amplitude distribution of CSD (a = 507 Hz, f = 204 Hz).

Z.M. Chen et al. / Applied Acoustics 87 (2015) 64–71 10

x 10 3.5 0.8 3

0.6 0.4

2.5

Y/m

0.2

2

0 1.5

-0.2 -0.4

1

-0.6

0.5

-0.8 -1

-0.5

0

0.5

71

motor, such hologram results do not truly reflect the radiation characteristics of the rotating or reciprocating machinery. Secondly, to select the CSD as the reconstruction variable, holography reconstruction was implemented at the cyclic frequencies corresponding to these frequencies, as shown in Figs. 15–20. The air compressor shield surface was chosen as the reconstruction plane. Figs. 15–20 show that the sound sources at 25 Hz, 17 Hz, 60 Hz, 237 Hz, 204 Hz and 93 Hz are mainly located at the foot of the driving motor, central part of the compressor, the air inlet filter of the compressor, the compression chamber of the compressor, the air outlet of the compressor and the inlet of the compressed air tank respectively. So the sound sources of the air compressor include the air inlet filter, the air outlet of the compressor and the inlet of the compressed air tank besides the driving motor and the compressor.

1

X/m Fig. 20. The amplitude distribution of CSD (a = 729 Hz, f = 93 Hz).

cannot be displayed on the hologram map. The reconstruction results based PNAH at a = 0, f = 188 Hz and d = ± 0.2 m are shown in Fig. 9. Obviously no energy distribution of the friction noise can be observed at the coupling. 4.2. Experiment on an air compressor The experiment was conducted on a CZ-20-C type air compressor. The rotating speeds of the driving motor and the compressor are 2990 r/min and 1000 r/min respectively. In order to attenuate sound reflection of the ground, a fabric carpet was laid on the ground. The measuring grid comprised 27(horizontal direction)  24(vertical direction) points with lattice spacing of 8 cm in both directions. The holographic plane was 10 cm away from the outer surface of the air compressor. The reference microphone was placed in front of the air compressor. The sampling frequency was 2048 Hz and the sampling time was 10 s. For eliminating the ground reflection, the carpet was laid on the ground. The experimental setup is shown in Fig. 10. The zero time delay CAF of sound pressure signal measured by the reference microphone is shown in Fig. 11 where the peak cyclic frequencies have been marked out. By analyzing the demodulation property of CAF, it is known that the modulating frequency and its harmonics appear in lower frequency band on the cyclic frequency axis, while double carrier frequency and sideband components of the modulating frequency appear in the higher frequency band. The compressor is the main source of the sound field radiated from the air compressor, and the rotating frequency of the compressor, 17 Hz, is also the main modulating frequency of the sound field. CSD slices at the peak cyclic frequencies 0 Hz, 17 Hz, 68 Hz, 219 Hz, 407 Hz, 507 Hz and 729 Hz were calculated respectively. The gathering slices of CSD are shown in Fig. 12. Fig. 12 indicates that the CSD is the power spectral density when a = 0. When a takes values other than zero, frequency components of 25 Hz, 17 Hz, 60 Hz, 237 Hz, 204 Hz, 93 Hz appear on the f frequency axis. There are sidebands with frequency spacing of 17 Hz surrounding each of these frequencies. First, to select the PSD as the reconstruction variable, a = 0, f = 17 Hz and a = 0, f = 204 Hz as the reconstruction frequency, air compressor shield surface as reconstruction plane. As shown in Figs. 13 and 14, if periodical time-variant properties were ignored, whether the compressor shaft frequency 17 Hz or the demodulation frequency 204 Hz as reconstruction frequency, the main noise sources of air compressor are located at the foot of the driving

5. Conclusions Second-order cyclic statistics have unique advantages in the analysis of cyclostationary signal. In this article, a PNAH technique for cyclostationary sound fields, known as CPNAH technique has been presented. Meanwhile, the gathering slices method for preprocessing of holographic data has also been presented. In the motor experiment, the cyclic frequency of 16 Hz produced by friction was found by analyzing cyclic autocorrelation function and CSD slice, and the distribution of friction acoustic noise was obtained by CPNAH reconstruction, which is impossible to be accessible by conventional PNAH technique. In the air compressor experiment, main sound sources of the air compressor, including the compressor, the driving motor, the air inlet filter, the air outlet of compressor and the inlet of the compressed air tank, were identified and localized successfully. The experimental results prove that both CPNAH technique and the gathering slices method are effective in analyzing cyclostationary sound fields. References [1] Randall RB, Antoni J, Chobsaard S. Comparison of cyclostationary and envelope analysis in the diagnostics of rolling element bearings. In: ICASSP, IEEE international conference on acoustics, speech and signal processing – proceedings, vol. 6; 2000. p. 3882–5. [2] Randall RB, Antoni J, Chobsaard S. The relationship between spectral correlation and envelope analysis in the diagnostics of bearing faults and other cyclostationary machine signals. Mech Syst Signal Process 2001;15(5):945–62. [3] Antoni J, Chobsaard S. Cyclostationary analysis of rolling-element bearing vibration signals. J Sound Vib 2001;248(5):829–45. [4] Antoni J, Bonnardot F, Raad A, Badaoui ME. Cyclostationary modeling of rotating machine vibration signals. Mech Syst Signal Process 2004;18:1285–314. [5] Gardner WA. Measurement of spectral correlation. IEEE Trans Acoust Speech Signal Process 1986;34:1111–23. [6] Brown WA. On the theory of cyclostationary signals. Ph. D dissertation, Department of Electrical Engineering and Computer Science, University of California, Davis; 1987. [7] Maynard JD, Williams EG, Lee Y. Near field acoustic holography I. Theory of generalized holography and the development of NAH. J Acoust Soc Am 1985;78(4):1395–413. [8] Veronesi WA, Maynard JD. Near-field acoustic holography (NAH) II: holographic reconstruction, algorithms and computer implementation. J Acoust Soc Am 1987;81(5):1307–22. [9] Wan Q, Jiang WK. Near field acoustic holography (NAH) theory for cyclostationary sound field and its application. J Sound Vib 2006;290:956–67. [10] Prakriya S. Blind identification of nonlinear systems based on higher order cyclic spectra. Ph. D. dissertation, University of Toronto (Canada); 1997. [11] Smith L. Blind channel identification and equalization using second-order cyclostationarity. Ph. D. dissertation, The Pennsylvania State University; 1996. [12] Jason FD. Time aliasing methods of spectrum estimation. USA: Faculty of Brigham Young University; April 2003. [13] Roberts RS, Brown WA, Loomis HH. Computationally efficient algorithms for cyclic spectral analysis. Signal Process 1991;4:38–49. [14] Timar PL, Fazekas A, Kiss J, Miklos A, Yang SJ. Noise and vibration of electrical machines. Amsterdam: Elsevier Science Publishers; 1989.