WO2015047065A1

WO2015047065A1 - Method and device for estimating the real and imaginary parts of the frequency response

Info

Publication number: WO2015047065A1
Application number: PCT/LV2013/000011
Authority: WO
Inventors: Vairis STRAUSS; Aldis KALPINS; Uldis LOMANOVSKIS
Original assignee: Latvijas Universitates Agentura "Latvijas Universitates Polimeru Mehanikas Instituts"
Priority date: 2013-09-27
Filing date: 2013-09-27
Publication date: 2015-04-02

Abstract

A method for estimating the real and imaginary parts of the frequency response for a physical system, the method comprising the steps of: acquiring values of the amplitude response of the system for a number of geometrically spaced frequencies; calculating the real part of the frequency response from the acquired values of the amplitude response by a first digital filtering; and calculating the imaginary part of the frequency response from the acquired values of the amplitude response by a second digital filtering. A device for carrying out the method is also proposed.

Description

Method and device for estimating the real and imaginary parts of the frequency response

Technical Filed The present invention relates to the field of estimating the real and imaginary parts of the frequency response of a physical system, which may be any object that could be described as a linear time- invariant or shift- invariant system having monotonic responses to impulse and step excitations.

Background Art Usually the real and imaginary part of the frequency response are determined by applying an excitation signal to the system, measuring the response signal of the system and then comparing the response signal to the excitation signal. Most up-to-date methods of determination or estimation of the real and imaginary part of the frequency response are based on performing the Fourier transform of the response and excitation signals, which is usually carried out by: correlation calculation, where the signals are multiplied by a complex exponential with integrating the product (U.S. Pat. No. 7,428,683 Sep. 23, 2008, U.S. Pat. No. 4,713,782 Dec. 15, 1987);

calculation by some version of discrete Fourier transforms (DFT) (U.S. Pat No.

4607216 Aug. 19, 1986, U.S. Pat. No. 4,991,128 Feb. 05, 1991).

These methods comprise the following operations performed according to the diagram given in Fig. 1 : a) generating an excitation signal to be applied to a system;

b) sampling the excitation signal and converting it by an Analog-to-Digital (A/D) converter into a digital signal;

c) sampling the response signal and converting it by an A/D converter into a digital signal; d) storing the both digital signals;

e) performing Fourier transform of the stored digital signals (by means of correlation calculation or with some versions of DFT);

f) dividing the Fourier transform of the response signal by the Fourier transform of the excitation signal.

Realization of all the mentioned signal digitisation and data acquisition procedures leads to high technical complexity of the structures supporting them, in particularly, in the cases, when frequency responses shall be determined over wide frequency ranges.

The procedures related to processing of the acquired data are complicated. Transforming the acquired time-domain signals into the frequency domain is a complicated digital processing, which requires powerful processors for fulfilling the Shannon sampling theorem [1] (sampling rate must be at least two time of the highest frequency of the signal). In the cases, when frequency responses shall be determined over wide frequency ranges, a large quantity of the acquired data must be processed. At the same time, since frequency responses are calculated at equally spaced points, the large redundancy of points appears in the high-frequency region. The necessity of performing these computationally burdensome calculations at digital processing of signals is a substantial drawback of the conventional methods for estimating the real and imaginary parts of the frequency responses. Thus the main problem is the complexity reduction in estimating the real and imaginary parts of the frequency responses.

Disclosure of Invention

The object of the invention is to reduce the complexity of estimating the real and imaginary parts of the frequency responses. According to the invention, an estimation of the real and imaginary parts of the frequency response is carried out through the amplitudes measurements of sinusoidal excitation and response signals at geometrically spaced frequencies by means of digital filtering.

In one aspect the invention provides a method for estimating the real and imaginary parts of the frequency response for a physical system, the method comprising the steps of: 2013/000011

3 a) acquiring values of the amplitude response of the system for a number of geometrically spaced frequencies;

b) calculating the real part of the frequency response from the acquired values of the amplitude response by a first digital filtering;

c) calculating the imaginary part of the frequency response from the acquired values of the amplitude response by a second digital filtering.

The step of acquiring values of the amplitude response may comprise generating sinusoidal excitation signals at geometrically spaced frequencies and applying the generated signals to the system; as well as measuring amplitudes of the sinusoidal excitation signals and measuring amplitudes of the response signals of the system; and storing the measured amplitudes of the sinusoidal excitation and response signals for determining values of the amplitude response of the system at geometrically spaced frequencies by dividing the measured amplitudes of the sinusoidal response signals by the measured amplitudes of the sinusoidal excitation signals. The main advantages of the proposed acquiring values of the amplitude response are as follows: it is easy to generate sinusoidal excitation signals at different frequencies;

geometrically spaced frequencies eliminate the redundancy of points of frequency response in the high-frequency region; if frequencies are spaced geometrically, a small number of sinusoidal excitation signals must be generated to cover frequency responses over wide frequency ranges; it is relatively easy to implement accurate measurement of amplitudes of the sinusoidal signals, for example, by rectifying a sinusoidal signal and measuring the rectified direct current (dc) signal; it is easy to obtain values of the amplitude response by dividing the amplitudes of the response signals by the amplitudes of the excitation signals, for example, if a sinusoidal excitation signal with the unit amplitude is used, the amplitude response is equal to the amplitude of the sinusoidal response signal.

Preferably, the digital filtering is carried out by convo luting the values of the amplitude response with the coefficients of a finite impulse response (FIR) filter.

FIR filters have the following advantages: they are computationally efficient algorithms working without employing numerical integration. They have the well-developed theory, see for example [1]. FIR filters may be constructed with the uniform structure and implementation in software and hardware and they can be modified very easy by changing filter coefficients without modification of the common structure or implementation of the algorithm in hardware and software. For correctly designed FIR filters no stability problems occur and they have the guaranteed performance, such as accuracy, sensitivity to noise, etc. Preferably, the frequencies are geometrically spaced with a common ratio of two.

For the common ratio of two, the Shannon sampling theorem is satisfied for geometrically sampled values of the amplitude responses (in the Mellin transform domain [2]); the common ratio of two is easy to be implemented by digital electronics; the common ratio of two results in small number of filter coefficients and small number of values of the amplitude response to be filtered, for example, 21 values of the amplitude response at frequencies spaced geometrically at the common ratio of two cover frequency range over more than 6 decimal decades. Compared with equally spaced samples, more than 10⁷ samples are necessary to cover this frequency range.

For the geometrically spaced frequencies with the common ratio of two, the digital filter with the following eight coefficients may be used for calculating the real part:

{-.00862702, 0.0474832, -0.143979, 0.480793, 0.759683, -0.185825, 0.0656697,

-0.0153284}.

For the geometrically spaced frequencies with the common ratio of two, the digital filter with the following eight coefficients may be used for calculating the imaginaiy part:

{0.109042, 0.0134620, 0.165268, 0.601015, -0.372957, -0.470905, 0.136201, -0.178838}.

Since a number of filter coefficients determines accuracy of the real and imaginary parts of the frequency response calculated by filtering, the larger number of coefficients, the higher is the accuracy. Eight coefficients provide relatively high accuracy with error of few percents. The method may be used for a variety of physical systems, such as mechanical, electrical, acoustical, etc. Network analysis, circuit simulation, material science, relaxation spectroscopy are just a few examples of the application areas for the method. In particularly the method, may be advantageously applied for estimation of complex dielectric permittivity of material. 1

5

In another aspect, the invention provides a device for performing the method, the device comprising: an excitation synthesiser for generating sinusoidal excitation signals at geometrically spaced frequencies to be applied to the system; a measuring unit for measuring amplitudes of the sinusoidal excitation signals and response signals of the system; a processing unit for calculating the real and imaginary parts of the frequency response of the system from the measured amplitudes of the sinusoidal excitation signals and response signals of the system.

The processing unit may comprise a first digital filter with the coefficients adapted to calculate the real part of the frequency response and a second digital filter with the coefficients adapted to calculate the imaginary part of the frequency response.

Alternatively, the processing unit may contain one digital filter with switching coefficients for calculating both real and imaginary parts in sequence.

Preferably, the digital filters have finite impulse responses.

An advantage of the device for performing the method is its simplicity, no sophisticated measuring appliances or powerful processors being required for its implementation.

In still another aspect, the invention provides a computer program adapted to carry out the method. That program may be executed on a suitable general-purpose computer or on a specialised computer.

The proposed invention reduces the complexity of the process of estimating the real and imaginary parts both at the data acquisition stage and at the data processing stage.

The complexity reduction in data acquisition is achieved by measuring a small number of amplitudes of sinusoidal excitation and response signals at geometrically spaced frequencies instead of digital recording of a large number of equally spaced samples of excitation and response signals in the time domain required by the sampling rate complying with the Shannon sampling theorem [1].

The complexity reduction in estimating the real and imaginary parts of the frequency responses is achieved by means of digitally filtering a small number of divisions of the measured amplitudes of the sinusoidal excitation and response signals instead of 1

6 performing the Fourier transform of a large number of samples of the excitation and response signals and dividing the appropriate Fou ier transforms.

Brief Description of Drawings

Fig. 1. Generalized block diagram of determining the frequency response according to the prior art methods.

Fig. 2. Generalized block diagram of a method for estimating the real and imaginary parts of the frequency response according to the invented method.

Fig. 3. Flowchart of the operations to calculate both the real and imaginary parts of the frequency response by filtering the amplitude response of the system

determined at geometrically spaced frequencies.

Fig. 4. Block diagram of a device, which uses the method proposed.

Fig. 5. Measuring circuit of a system for measurement of complex dielectric permittivity of a material.

Fig. 6. Block diagram of a measuring unit.

Fig. 7. Processing unit of Fig. 4 in more detail.

Fig. 8. Diffuse magnitude response of the ideal filter for determining the real part of the frequency response (shaded area) and its approximation by the magnitude response of digital filter with 8 coefficients at the common ratio of two (dashed line).

Fig. 9. Diffuse magnitude response of the ideal filter for determining the imaginary part of the frequency response (shaded area) and its approximation by the magnitude response of digital filter with 8 coefficients at the common ratio of two (dashed line).

Detailed Description of Invention

Fig. 2 describes a generalized architecture of the proposed method for calculating the real and imaginary parts of the frequency response of the system. As compared to the Fig. 1 , there are substantial differences in the approach to both the A/D conversion and data processing. V2013/000011

The invented method for estimating the real and imaginary part of the frequency response of the system comprises: a) generating sinusoidal excitation signals at geometrically spaced frequencies to be applied to the system;

b) measuring amplitudes of the sinusoidal excitation signals;

c) measuring amplitudes of the response signals of the system;

d) storing the measured amplitudes of the sinusoidal excitation and response

signals;

e) determining values of the amplitude response of the system at geometrically spaced frequencies by dividing the measured amplitudes of the sinusoidal response signals by the measured amplitudes of the sinusoidal excitation signals;

f) calculating the real part of the frequency response of the system by digital filtering the values of the amplitude response;

g) calculating the imaginary part of the frequency response of the system by

digital filtering the values of the amplitude response.

The data acquisition process represented by the steps a) to e) above may be carried out separately from the calculating steps and may be carried out by any other methods known in the art.

The calculation of both the real and imaginary parts of the frequency response is composed of the following sub-operations (Fig. 3 and Fig 7):

(1) Fetching the first (N- l) signal samples from Data Input Buffer and entering into the Shift Registers, where N is the number of coefficients of digital filter.

(2) Fetching the next signal sample from Data Input Buffer and shifting the content of Shift Registers by one sample.

(3) Multiplying the sample value of output of each stage of Shift Registers by the corresponding coefficient of Digital Filter.

(4) Summing the digital products of signal samples and coefficients.

(5) Storing the value of filter digital output signal into Data Output Buffer, and initiating the next cycle of filtering. 1

8

In case of measuring complex dielectric permittivity of a material, the system may represent a parallel plate capacitor with inserted material under test. The measuring circuit of such system is shown in Fig. 5.

The parallel plate capacitor 1 produces excitation of the form of electrical field to the material under test 2. The response signal is obtained as voltage drop on measuring resistance 3.

The Processing unit carries out the digital filtering of the magnitude response according to the proposed method. All the other blocks shown in Fig. 4 are involved to acquire the values of the amplitude response sampled according to geometrical progression in the frequency domain.

The Excitation synthesizer generates sinusoidal excitation signals with frequencies in a geometric progression. Excitation signals are applied to the system to excite response signals.

Both excitation and response signals are applied to a multiplexer Ml, which is used to select either signal for further processing by Measuring unit. The digitised values of measured amplitudes are then passed to Processing unit for storage and processing. The values of the amplitude response of the system are computed by dividing the measured amplitude values of the response signal by the appropriate measured amplitude values of the excitation signal. The real and imaginary parts of the frequency response of the system are calculated by digital filtering the values of the amplitude response of the system.

Fig. 6 shows a generalized block diagram of the Measuring unit.

The Measuring unit measures amplitudes of sinusoidal excitation and response signals. It contains Variable gain amplifier, Rectifier, Integrator, and A/D Converter.

In order to make full use of the dynamic conversion range of A/D Converter, a sinusoidal signal at first goes through a digitally controlled Variable gain amplifier. Rectifier converts the electrical sinusoidal signal to some form of pulsating dc signal. Afterwards those dc pulses are integrated with respect to time producing a dc signal, which represents the amplitude of sinusoidal signal. The resulting dc signal is measured with an A/D converter for further processing by the Processing unit. V2013/000011

9

Fig. 7 shows the hardware implementation of the Processing unit.

The Processing Unit includes Data Input Buffer, Divide Unit, Digital filter for calculating real part, Digital filter for calculating imaginary part, and Data Output Buffer. The both Digital filters have the same architecture and differ by the content of Coefficient storage. Input signal samples for the Processing Unit are digitised values of the measured amplitudes, which are loaded into Data Input Buffer. To obtain the values of the amplitude response of the system, the content of Data Input Buffer is normalized by Divide unit, where the values of amplitudes of the response signals are divided by the appropriate values of amplitudes of the excitation signals. Initially the first (N- 1 ) amplitude response values are loaded into Shift Registers of the both Digital filters. Then the next amplitude response value is fetched from Data Input Buffer and the content of Shift Registers is shifted by one digital sample. Thereafter the digital sample value of output of each stage of Shift Registers is multiplied by the corresponding coefficient of Digital filter. The value of Digital filter output is obtained by summing the products of Digital signal samples and coefficients. To obtain the next output value of Digital filter the next amplitude response value is fetched from Data Input Buffer and the content of Shift Registers is shifted by one digital sample as well as the mentioned multiplying and summing operations are repeated. The output values of Digital filters are stored into Data Output Buffer.

Thus the proposed method makes it possible to avoid bulky digital recording samples of the excitation and response signals dictated by the Shannon sampling theorem [1], and significant simplification of the digital processing by a much simpler operation of digital filtering by filters having short impulse responses instead of performing Fourier transformations.

Determination of the filter parameters

In general, the frequency response J{d) of physical systems is a complex function J( o) = J'(o) + jJ"((o) , where J'{(o) is the real part, J"(a>) is the imaginary part, j = V- . The modulus of the frequency response

is named amplitude response. The real part J'(ct>) of physical systems with monotonic responses varies monotonically from value J₀ at zero frequency to value J_∞ at infinity frequency as the frequency increases, while the imaginary part J"{co) has zero value at zero and infinite frequencies and passes through a maximum(s) at some frequency(ies) between zero and infinity.

The real and imaginary parts of the frequency response of causal physical systems, such as materials, electrical and mechanical devices, electrical networks, etc., are not wholly independent but are linked by a special form of Hilbert transforms, which are termed the Kramers-Kronig (KK) relations [3] allowing to calculate J"(co) from J'(co) , and vice versa.

Amplitude response of the system | J(co) | according to (1) also contains information about the real and imaginary parts, thus, it would be potentially possible to calculate the real and imaginary parts from the amplitude response.

In the present invention, the real and imaginary parts are calculated by digital filtering from the values of the amplitude response (1) determined at geometrically spaced frequencies. In the case when the imaginary part is small compared to the real part, the amplitude response of the system approaches to the real part

\ J(a>) |→J» , and the imaginary part can be calculated from the amplitude response by the KK relation [3] 1 1

while the real part is approximately equal to the amplitude response

J'((o) \ J(a)) \ . (3)

As shown in [4], KK relation (2) may be interpreted as a linear shift-invariant system or an ideal filter operating on a logarithmic frequency scale having the following frequency response in the Mellin transform domain

where j = V-T , and ]μ is the Mellin transform parameter.

In its turn, determination of the real part from the amplitude response according to (3) can be also treated as a filtering problem with an all-pass filter having a constant unity frequency response

The phenomeno logical theories [5,6] state the upper bound for the imaginary part of the frequency response of an elementary causal system as A/COT . .. corresponding to the following real part

J_m' co) = J_{∞ +} -^ , (5)

I + CO % where AJ = J_Q - J_∞

It can be shown that determination of the upper bound or maximum imaginary part (4) from the amplitude response (1) can be also treated as a filtering problem with a filter having the frequency response Η' (μ) = _ V± M (6) m^ax 2^^μ Γ{2 + βμ)ΓΗμ∞*απμΙ2Υ where Γ(.) is gamma function.

Similarly, determination of the real part (5) from the amplitude response (1) can be treated as a filtering problem with a filter having the frequency response _H, _{(u) =} πΓ{2 + ]μ)

"^m 2¹-"¹ Γ (2 + ]2μ)Γ (- /μ) ΰη(-]πμ 12) -

In Fig. 8, amplitude responses | H_m' _in (μ) | and | H_m' _ax (μ) | are shown for the filters for determining the limiting real parts (3) and (7) from the amplitude response (1). For determination of the real parts from the amplitude response (1) in all other cases, the amplitude responses of the filters shall lie in the area (shaded) between | H_m' _ia (μ) | and I H_m' _ax (μ) I . Therefore, the problem of determination of the real part from the amplitude response of the system may be interpreted as a filtering task with a diffuse amplitude response | Η'(μ) | bounded by the limiting amplitude responses | H_m' _in (μ) | and | H^ (μ) | .

Similarly, in Fig. 9, the amplitude responses | H",,,, (μ) | and j H_m"_M (μ) | are shown for the filters for determination of the limiting imaginary parts (2) and (4). In the same manner as for the real part, the problem of determination of the imaginary part J"(co) from the amplitude response (1) may be interpreted as a filtering task with a filter with diffuse amplitude response | Η μ) \ bounded by the limiting amplitude responses | H^_n (μ) | and

Digital filters for calculating the real part J'(co) and the imaginary part J"((o) from the amplitude response | J( o) | may be implemented in the form with finite impulse responses [1]

(tf-l)/2

'(¾) = ∑ h[n]x{co_m - ω„) , (8) 2013/000011

13 where y( _m) is output function representing the real and imaginary parts of the frequency response (i.e. y{co_m) = J'(fi>,„) and y(co_m ) = J"(co_m) ), *(&>„) is the input function representing the magnitude response ( x(o) = J(ft>„) ), and N is a number of filter coefficients. Since uniformly sampled data on the logarithmic scale manifest as the samples distributed according to a geometric progression in the linear scale [4], frequencies <¾ for the digital filter (8) are geometrically distributed io„= co₀ ", « = 0,± 1,± 2...., (9) where q > 1 is common ratio. Therefore the digital filter (8) executes an algorithm [4]:

(W-l) / 2

y(<o₀q^m) = ∑h[n]x( _oq"-") (10) and has the periodic, with the period of 2π I In q , frequency response in the Mellin transform domain

H(0 =∑h[n]zxp(- !^ \_n q) . (1 1)

« Despite that ideal amplitude responses | Η μ) | and | Η"(μ) | are diffused, they can be approximated by single discrete impulse responses h[n] . The filter coefficients, in principle, can be obtained by conventional design methods of digital filters [1]

implementing inverse discrete Fourier transform (IDFT) of some chosen frequency responses, which are bounded by the limiting frequency responses | H^,, (μ) | and

I HI, (μ) I for calculating the real part J c ) and | H"_u„ (μ) \ and | H^ (μ) | for calculating the imaginary part. However, the best results give the method of designing digital filters based on the system identification principle [4, 7], where the filter coefficients h[n] are determined from the equation (8) or (10) by using chosen pair(s) of the exact input and output functions x(co) and y(co) . In Table I below, exemplary coefficients are given for filters for calculating the real and imaginary parts containing N - 8 coefficients and operating at a common ratio of q - 2. In Fig. 7, the filter coefficients h'[n] and /?"[«] are used re-indexed to run from 1 to 8. The appropriate amplitude responses (11) are shown in Fig. 8 and Fig. 9 in dashed lines. Table I

Industrial Applicability

The proposed method and device may be applied for estimating the real and imaginary parts of the frequency responses of any physical systems, which have monotonic responses to impulse and step excitations, including but not limited to, materials, electrical and mechanical devices, electrical networks, etc.

References

Patent documents 1. U.S. Pat. No. 7,428,683 Sep. 23, 2008,

2. U.S. Pat. No. 4,713,782 Dec. 15, 1987 LV2013/000011

15

3. U.S. Pat. No. 4,607,216 Aug. 19, 1986

4. U.S. Pat. No. 4,991,128 Feb. 05, 1991

Non-patent documents

[1] Oppenheim, A. V., Schafer, R. V. Discrete-Time Signal Processing; 2nd. Ed.;

Prentice-Hall International, New Jersey, 1999.

[2] Bertrand, J., Bertrand, P., Ovarlez, J.-P. Chapter 11. The Mellin Transform. In: The

Transform and Applications Handbook, 2. ed. Ed. A. D. Poularikas, CRC Press, Boca

Raton, 2000.

[3] Nussenzveig, H. M. Causality and Dispersion Relations; Academic Press, NY, 1972.

[4] Shtrauss, V. FIR Kramers-Kronig transformers for relaxation data conversion, Signal

Processing 2006, 86, 2887-2900.

[5] McCrum, N. G., Read B. E., Williams G., Anelastic and Dielectric Effects in

Polymer Solids, J. Wiley and Sons, London, 1967.

[6] Tschoegl, N. W. The Phenomenological Theory of Linear Viscoelastic Behavior,

Springer- Verlag, Berlin, 1989.

[7] Shtrauss, V. Functional conversion of signals in the study of relaxation phenomena,

Signal Processing, 1995, 45, 293-312.

Claims

1. A method for estimating the real and imaginary parts of the frequency response for a physical system, the method comprising the steps of:

(a) acquiring values of the amplitude response of the system for a number of geometrically spaced frequencies;

(b) calculating the real part of the frequency response from the acquired values of the amplitude response by a first digital filtering;

(c) calculating the imaginary part of the frequency response from the acquired values of the amplitude response by a second digital filtering.

2. The method of claim 1 wherein the step (a) comprises generating sinusoidal excitation signals at geometrically spaced frequencies and applying the signals to the system.

3. The method of claim 2 wherein the step (a) further comprises measuring amplitudes of the sinusoidal excitation signals and measuring amplitudes of the response signals of the system.

4. The method of claim 3 wherein the step (a) further comprises determining values of the amplitude response of the system by dividing the measured amplitudes of the response signals by the measured amplitudes of the sinusoidal excitation signals.

5. The method of claim 4 wherein the step (a) further comprises storing the measured amplitudes of the sinusoidal excitation and response signals.

6. The method of any of claims 1-5 wherein the digital filtering is carried out by convoluting the values of the amplitude response with the coefficients of a finite impulse response filter.

7. The method of claim 6 wherein the frequencies are geometrically spaced with a common ratio of two.

8. The method of claim 7 wherein in the step (b) the digital filter with the following eight coefficients is used: {-0.00862702, 0.0474832, -0.143979, 0.480793, 0.759683, -0.185825, 0.0656697, -0.0153284}.

9. The method of claim 7 or 8 wherein in the step (c) the digital filter with the following eight coefficients is used: {0.109042, 0.0134620, 0.165268, 0.601015, -0.372957,

-0.470905, 0.136201, -0.178838}.

10. The method of claim 1-9 used for estimating the complex dielectric permittivity of a material.

11. A device for estimating the frequency response for a physical system comprising: an excitation synthesiser for generating sinusoidal excitation signals at

geometrically spaced frequencies to be applied to the system; a measuring unit for measuring amplitudes of the sinusoidal excitation signals and response signals of the system; a processing unit for calculating the real and imaginary parts of the frequency response of the system from the measured amplitudes of the sinusoidal excitation signals and response signals of the system.

12. The device according to claim 11 wherein the processing unit comprises a first digital filter adapted to calculate the real part of the frequency response and a second digital filter adapted to calculate the imaginary part of the frequency response.

13. The device according to claim 12 wherein the digital filters have finite impulse responses.

14. A computer program adapted to carry out the method according to any one of the claims 1-9.