Digital Frequency Discriminator, for Extracting a Sinusoidal Signal

The present invention relates to a digital frequency discriminator intended for extracting from a carrier at a frequency of f.sub.o at least one low frequency sinusoidal signal, of frequency F, frequency modulating the said carrier with a modulating frequency excursion .DELTA.F, the values of f.sub.o and .DELTA.F being known.

The present invention further relates to the use of a discriminator for determining the phase known as the reference phase contained in the composite VOR signal, received by an on-board VOR-ILS equipment, by point by point restoration of the 30 Hz sinusoidal signal giving the reference phase.

Most of the frequency discriminators used up to the present time are analog discriminators. They are used in particular in radio navigation. For example, in a radio altimeter of the linear frequency modulated continuous wave type, it is a known technique to hold the subtractive beat frequency f.sub.b between the transmitted wave and the received reflected wave at a constant set frequency f.sub.bo of predetermined value by including such a discriminator in the feedback system which controls the slope of the modulation line of the transmitted signal. This discriminator which has the function of providing, according to a substantially linear characteristic, the difference between the frequency f.sub.b and the set frequency, is basically formed of two processing systems each of which receives the signal to be discriminated. Each processing system includes a band pass filter, centered on f.sub.bo -.delta.f (and centered on f.sub.bo +.delta.f) respectively), followed by an absolute value detector. The two signals thus obtained are subtracted from each other and the resultant signal is filtered through a low pass filter, An S-shaped characteristic amplitude-frequency curve is thus obtained whose central section, on either side of a phase shift equal to .pi./2, is substantially linear. From the work "Digital Signal Processing" by V. Capellini and A. G. Constantini, published by Elsevier Science Publishers B.V. (North Holland 1984), pages 632 to 637, the technique is known of digitally transposing the structure of the analog discriminator described above. The six components which form this structure are thus produced digitally, which obtains certain advantages such as ease of adjustment during mass production for example. In this way the previously mentioned S-shaped characteristic curve is again obtained. This digital discriminator has disadvantages: its linearity is limited. Also, the filtering thus achieved requires a large computing capacity and the result obtained, even if it is acceptable for a continuous frequency spectrum having a leading edge which it is sought to discriminate as is the case for the beat signal of a radio altimeter, has a lower performance when the discriminator of a single frequency line is desired, for example by demodulation of a VOR-ILS signal.

A purpose of the invention is to produce a digital discriminator enabling the reconstruction, with the help of a digital process, of at least one low frequency signal transmitted on a frequency modulated carrier.

Another purpose of the invention is to produce a digital discriminator having excellent linearity.

Yet another purpose of the invention is to reduce the amount of computations carried out by the digital discriminator.

These purposes are achieved due to the fact that the digital discriminator of the type specified in the opening paragraph according to the invention is characterized in that it receives, with or without prefiltering, digital samples x.sub.n of the said carrier modulated at a frequency f.sub.e such that f.sub.e =4 f.sub.o, and in that it includes

means of accumulation for accumulating between two successive dates t.sub.q and t.sub.q+1 a number N, which can be set, of successive digital samples x.sub.n which form an analysis window. First means of storage for storing coefficients R.sub.01, R.sub.02, R.sub.11, R.sub.12 and R.sub.22 defined by: ##EQU3##

means of formatting the said coefficients, first means of computing the quantity a.sub.1 defined by: ##EQU4## the quantity a.sub.1 being similar, to the nearest multiplicative coefficient, to the ordinate of a point on the curve representing the said sought low frequency sinusoidal signal, and in that the value of a.sub.1 is transferred into second means of storage, after which the next point of the said low frequency sinuosoidal signal is determined from the N samples which follow the data t.sub.q+1.

The computations carried out in order to perform a point by point restoration of the said low frequency sinusoidal signal result from a choice of a simple second order mathematical model and simplifying assumptions explained in detail hereafter. The purpose of formatting the coefficients is to examine, for each restored point, whether the coefficients R.sub.01, R.sub.02, R.sub.11, R.sub.12 and R.sub.22 have a sufficiently large value for the determination of the considered point to be able to be declared valid and, in the affirmative, to facilitate the subsequent computations carried out on these coefficients for determining the quantity a.sub.1.

A preferred embodiment of the discriminator according to the invention is characterized in that it also includes second means of computing the quantity a.sub.2 defined by: ##EQU5## the quantity a.sub.2 being used as a convergence indicator and, in that for this purpose the said digital frequency discriminator includes means of comparison of the quantity a.sub.2 with two thresholds s.sub.1 and s.sub.2 whose values can be set such that the quantity a.sub.1 is only validated and transferred into the said second means of storage if: s.sub.1 <a.sub.2 <s.sub.2.

The value to be given to s.sub.1 and s.sub.2 depends on the amplitude of the useful signal to be processed. In the moist frequency case, it is a microwave signal containing noise and, below a certain value of the signal/noise ratio, it is considered that the discrimination of a low frequency signal is no longer possible. The determination of the thresholds s.sub.1 and s.sub.2 therefore results from the test of the equipment depending on the power of the received useful signal and the power of the noise mixed with this signal.

A particularly useful use of the discriminator to a VOR signal as mentioned in paragraph 2 of the present document is according to the invention characterized by a sampling frequency, equal to 13280 Hz, applied to a carrier of 3320 Hz obtained by filtering the 9960 Hz carrier of the VOR signal through an analog one-way filter of order less than or equal to three whose pass band is between 0 and 10600 Hz and whose transition band extends between 10600 Hz and 15960 Hz, by sampling the signal obtained at the output of the one-way filter at a frequency of 26260 Hz and by filtering the samples through a half-band filter with 4 distinct coefficients whose high pass output is connected to the input of the said digital discriminator.

By thus lowering the sampling frequency by reducing the frequency of the carrier to one third of its value, the volume of computation that must be carried out by the digital discriminator is consequently lowered.

The following description given with reference to the appended drawing figures, all given by way of example, will give a better understanding of how the invention can be carried out.

FIG. 1 is a block diagram of the signal processing system which includes a digital discriminator according to the invention.

FIG. 2 shows curves of estimation variance as a function of the signal/noise ratio of the output samples of the discriminator for various analysis windows and depending on whether or not prefiltering is included.

FIG. 3 shows the spectrum of a composite VOR signal.

FIGS. 4 and 5 each show a processing system that can be used for an analog VOR signal up to the production of the input samples e.sub.n of FIG. 1.

FIG. 6 shows a band separating digital filter preferably used in the processing system of FIG. 5.

FIG. 7 is a general flowchart explaining the functioning of the discriminator according to the invention.

FIGS. 8a and 8b form a detailed flowchart explaining the functioning of the discriminator according to the invention.

FIG. 9 is the computing flowchart of an element of the series v.sub.q associated with the series v.sub.q.

The exact technical problem to be resolved consists in finding, using a digital process, a low frequency signal y(t) of frequency F transmitted by frequency modulation on a carrier at the frequency f.sub.o. The low frequency modulating signal is written in simplified form as follows:

in which

a is the amplitude

.omega.=2.pi.F is the pulsation

.phi. is the phase at the origin.

In discrete form at the time of sampling i, the formula (1) can be written:

f.sub.e being the sampling frequency, the inverse of the sampling period .DELTA.T.

The frequency modulation carried out on the carrier with a frequency excursion .DELTA.F by y(t) is characterized as follows: the instantaneous frequency of the modulated carrier f.sub.in is expressed:

or again, in digital form:

The phase at time t, .PHI.(t) is therefore:

which can be written in discrete form at the sampling time n: ##EQU6## It will be noted that the phase at time n is connected to that at time (n-1) by a simple expression:

The expression of the frequency modulated carrier is therefore:

or, if this signal is coded in digital form: ##EQU7## where b is the amplitude, and .alpha. is the phase at the origin.

It is desired to restore the successive y.sub.n i.e. the modulation signal. The basic assumption of the digital processing process of the signal explained below consists in making the approximation that the instantaneous frequency f.sub.n (see formula 2) does not vary over a short time interval with respect to: T=1/F. In other words, the frequency f.sub.n for n such that: N.sub.o -N/2.ltoreq.n.ltoreq.N.sub.o +N/2 is, to a first approximation, equal to a constant f if N is small. Under these conditions, the problem raised is summarized by the point by point estimation of a sinusoid. Therefore, if f is the value of frequency to be estimated and .OMEGA.=2.pi.f.DELTA.T is the associated pseudo-pulsation; the following can be written to a first approximation:

because of replacing ##EQU8## By virtue of the trigonometric formula: ##EQU9## where p=n.OMEGA.+.alpha. and q=(n-2).OMEGA.+.alpha., the following emerges: sin (n.OMEGA.+.alpha.)+sin{(n-2).OMEGA.+.alpha.}=2 cos .OMEGA. sin{(n-1).OMEGA.+.alpha.} or again

The signals are received in a noisy environment and the problem consists in determining .OMEGA., and subsequently y.sub.n by virtue of the formula (2), under extremely variable conditions of signal to noise ratio S/N, this S/N ratio varying from 50 dB to 0 dB for example, from simple observations. For this, it is considered that each sample has noise and that at time n.DELTA.T:

is therefore obtained, x.sub.n being the digital sample of the input signal of the discriminator and b.sub.n the quantity of noise associated with this sample; this gives:

The expression (5) is then written:

or more precisely:

expressions in which p is the modulus of a complex number, close to unity, and .epsilon..sub.n represents the quantity of noise coming from the original noise and filtered through a band pass filter centered on f:

The recurrence, relationship between three successive samples (6) or (6)', represents a second order modelling of the system which corresponds to the detection of a single modulation line in the input signal of the discriminator. It will be noted that, for an input signal with several modulation lines, it is also possible, by means of a more complex modelling of the system, by an order greater than 2, to achieve the discrimination of these multiple lines. However, the following description will keep to the discrimination of a single line.

In order to estimate .OMEGA. the use of a linear prediction expression is proposed, of the type:

the expression (8) being able to be identified term by term, with the expression (6) or (6)'.

The noise .epsilon..sub.n is considered white, gaussian and therefore not correlated with the sought useful signal y(t). The estimation depends on the accumulation of the quantities .epsilon..sup.2.sub.n over a certain analysis window of N successive samples, which represent each time a point in the sought sinusoid, and on the minimization of the mean square error thus obtained, that is: ##EQU10##

The minimization of E is expressed as: ##EQU11## By development, two equations to be solved simultaneously are therefore obtained: ##EQU12## By putting ##EQU13## the above set of equations (11) is more simply expressed:

In the present case, the quantity a.sub.2 must equal 1 and p.sup.2 respectively if the expressions (6) and (8) or (6)' and (8) respectively are identified, but it is preferred to ignore this condition on first analysis and to use a.sub.2 as what may be called a convergence indicator. This means that if the computed value of a.sub.2 is too different from 1, it can be considered that there is too great a divergence and that the corresponding determination of y.sub.n must be rejected as being incorrect. Finally, this gives: ##EQU14##

The determination of the quantities a.sub.1 and a.sub.2 is subject to particular operating conditions and certain characteristics are favourable to such an estimation:

Firstly, the coefficients R.sub.11 and R.sub.12 represent the autocorrelation coefficient R.sub.0 which represents the power of the received signal and can be estimated over a wider horizon than the actual analysis window N on condition that there is standardization. For example, if N=64 is chosen, R.sub.0 can be estimated over a window of 128 or of 256 points. The estimation window of R.sub.0 must, however, remain rather narrow so that, for the time corresponding to this window, it can be considered that the signal has not been subject to a fluctuation in power.

Also essentially for reasons of simplicity of computation, a sampling frequency such that: f.sub.e =4f.sub.0 is preferably chosen. This choice also leads to optimum accuracy as demonstrated by the autoregressive methods in linear prediction theory. Thirdly, applications for which the frequency excursion of the modulating frequency .DELTA.F is small compared with the carrier frequency f.sub.o are desired, for example by a value of .DELTA.F of the order of 500 Hz for a value of f.sub.o of the order of 10,000 Hz. Finally the pass band of the modulating signal is considered as being small with respect to the carrier frequency f.sub.o.

Under these conditions, the digital processing to be carried out consists in dividing the observation signal, that is to say the modulated carrier, into adjacent sections each containing N samples. Each section gives rise to the determination of a point of the sought sinusoid. For each section of the signal, the approximation is made that the frequency of the observation signal remains constant and: ##EQU15## is obtained.

For f.sub.e =4f.sub.o and for R.sub.0 standardized at 1, it can be verified over the trigonometric circle that:

R.sub.01 tends to 0

R.sub.12 tends to 0

and

R.sub.02 tends to -1.

From this it is deduced that: a.sub.2 .apprxeq.1, which must be verified by measurements and computations. Similarly, ##EQU16## is obtained, and a.sub.1 tends to 0, the value of a.sub.1, although small, not however being negligible. If reference is made to the formulae (6) and (8),

may be written.

The fact that cos .OMEGA. is close to 0 can be used in a useful way by considering that .OMEGA. tends to .pi./2 and therefore: ##EQU17## form formulae (16), (17) and (18) it results that: ##EQU18## Now, the coming together of the formulae (3) and (4) under the assumption that: f.sub.e =4f.sub.o enables the following to be written: ##EQU19##

If, in this last formula, .OMEGA. is replaced by its value given by formula (19) ##EQU20## is obtained, in which formulae y.sub.n is the sought value of the modulating signal for the section concerned. Starting from formulae (6)' and (8), this results in a slightly more complicated but more rigorous computation at: ##EQU21##

In practice it is preferable to use expression (22)' rather than expression (22) particularly as the computed value of a.sub.2 moves further away from the value 1. Therefore, firstly the value of a.sub.2 is calculated, using formula (15) and only if this value is sufficiently close to 1, which indicates correct convergence, and then the values of a.sub.1 and y.sub.n are in their turn calculated using the formulae (16) and (22) or (22)' respectively.

It should be noted that insofar as the frequency excursion .DELTA.F is no longer negligible with respect to the carrier frequency f.sub.0, cos .OMEGA. no longer tends to 0 and consequently it becomes preferable to proceed with the extraction of .OMEGA. using the arc-cos function.

It is also posssible to compute the noise power from the formulae (8) and (9); this gives: ##EQU22## where ##EQU23##

According to the set of equation 12, ##EQU24## is obtained, .sigma..sup.2 being the noise power which can be estimated by means of expression (24).

As regards the embodiment, it will be noted that the expression (22) can be applied to an observation section of N points by computing the quantities R.sub.01, R.sub.02 and R.sub.12 using the expressions: ##EQU25## R.sub.0 being computed over a horizon which can be larger than N, namely L.

Another possibility consists in making a point by point recursive estimation using:

This discriminator according to the invention is designed, preferably, to function with signal to noise ratio S/N conditions that are widely variable, for example between 0 and 50 dB. When the S/N ratio drops below 15 dB, the functioning of the system thus described can be further improved a little by carrying out a prefiltering about the useful discrimination band, i.e. about the modulated carrier before discrimination about the frequency f.sub.o, on the assumption that: f.sub.o =f.sub.e /4. FIG. 1 shows the processing system of the signal including the digital discriminator according to the invention. The actual discriminator is preferably produced by means of a microprocessor as described below for the application of the invention to a composition VOR signal, such as shown in FIG. 3. The frequency modulated signal z(t) is supplied to a sampler-blocker 1 of a known type which samples at the frequency f.sub.e. The samples are then supplied to an analog/digital converter 2 whose output is the source of digital samples in serial form e.sub.n. A two-position switch 3, that can be controlled by a control signal, on a conductor 4, depending on the value of the S/N ratio (this signal having a different logic value above and below a threshold S/N=S'), enables the samples e.sub.n to be transmitted in the form x.sub.n to the input of the digital frequency discriminator 5 either directly or through a filter 6 which is a band pass filter on the frequency f.sub.o. An example of elementary band pass filtering centered on f.sub.o =f.sub.e /4 is given by the following expression:

This filter, although rudimentary, enables the noise power to be attenuated by means of simple operations such as addition and subtraction. The estimation variance of the modulation signal is thus improved for low values of the S/N ratio. The table 1 below gives an appreciation of the estimation standard deviation results for a modulating signal of amplitude a=1 in the case of observation (analysis) windows of 24, 48 and 96 points, for various values of the S/N ratio expressed in steps of 10 dB and with (A) or without (S) pre-filtering for each of these values.

FIG. 2 is a graph of the estimation standard deviation ESD given in the table, which represents the estimation variance as a function of the S/N ratio, the curves in full line indicating the discrimination without prefiltering and the curves in broken line indicating the discrimination with prefiltering; the number N of points in the observation window is indicated in brackets. This figure shows the effectiveness of the band pass filter between 15 dB and 0 dB where the standard deviation is particularly improved. Within this latter range of S/N ratio, the gain obtained on the S/N ratio with such a filter can be evaluated at 10 dB. FIG. 2 also shows that the number N must not be too small either as there is a risk of having too large an estimation standard deviation for the lowest values of the S/N ratio. Given the chosen values for various frequencies indicated above it is noted that the value N=24 would be too small while the values N=48 or N=96 are suitable. It is also possible to choose the value N=64 which, because it is a power of two, has the advantage of simplifying the computations.

It is proposed to apply the principle described above to the search for the reference phase in the VOR signal received by an on-board VOR-ILS equipment. The spectrum of the received composite signal A(FREQ.) is shown in FIG. 3. It contains a component at 30 Hz giving the variable phase, an audio band from 300 to 3400 Hz (or an identification signal at 1020 Hz), and a 9960 Hz carrier frequency modulated with an excursion of 480 Hz by the 30 Hz signal giving the reference phase. It is this frequency modulating 30 Hz signal that it is sought to restore.

A first solution for the application of the principle described above consists in sampling the VOR signal at a frequency of: 4.times.9960=39840 Hz, and then in separating the low sections of the spectrum (30 Hz component and audio signal) and the high sections (frequency modulated signal) and finally in applying the discrimination procedure to the 9960 Hz carrier. This latter frequency is equal to f.sub.e /4 and the computation procedure indicated above correctly applies. The processing system which converts the analog VOR signal into the digital signal e.sub.n is shown in FIG. 4. It is formed by the succession of a sampler-blocker 1, an analog-digital converter 2, and a digital separator filter 7, of known type, whose high pass output, referenced PH, is the source of the signal e.sub.n, at the same frequency as the input signal of the band separator 7. For the chosen frequency of f.sub.e at 39840 Hz, the volume of computation is large; it is advantageous to lower the sampling frequency.

Consequently, a second solution consists in starting by choosing a sampling frequency f.sub.e1 such that:

Before the sampling operation, the signal must in this case pass through an analog one-way filter in order to limit the maximum frequency of the spectrum to a value below 13280 Hz (cut-off frequency: f.sub.c =f.sub.e1 /2). As the high part of the VOR signal is frequency modulated with a nominal frequency excursion of 480 Hz and because of this extends to 10600 Hz, this one-way filter must have a pass band of 0 to 10600 Hz. The transition band is: 2(f.sub.e1 /2-10600)=5360 Hz and the attenuated band therefore extends from 10600+5360=15960 Hz to infinity. The order of the filter is determined by the ratio between the transition band and the cut-off frequency, namely: 5360/13280, giving an order of less than or equal to three. The operation which follows the sampling consists in separating the low and high sections of the spectrum. To do this, a digital band separator filter which is preferably a half band filter with 4 separate coefficients such as is shown in FIG. 6 is used. The processing system for the signal ahead of the digital discriminator is then the one shown in FIG. 5, consisting of a low pass analog one-way filter 8, a sampler-blocker 1, an analog/digital converter 2 and a digital band separator filter 9 shown in more detail in FIG. 6. The filter 9 is a finite pulse response halfband transverse filter, of known type, it operates alternately on even samples k.sub.n (n=2p) and odd samples k.sub.n (n=2p+1), which is represented in FIG. 6 by a switch. This structure enables it to carry out a subsampling by two. It receives samples k.sub.n at the rate of 26560 Hz and supplies at its high pass output PH samples e.sub.n at the rate of 13280 Hz. The structure of the filter in FIG. 6 is designed to perform the following equation: ##EQU26##

The delays z.sup.-1 shown in FIG. 6 are at the output rate, that is f.sub.e1 /2.

The response is therefore symmetrical on either side of a main coefficient and equal to 1 with the odd coefficients only, the even coefficients being zero. Because of the sub-sampling by two, the PH output is transposed by 13280 Hz, that is f.sub.e1 /2, and the new carrier before discrimination is at a frequency of 3320 Hz (13280-9960=3320 Hz); in effect, the sub-sampling by two gives, for the spectrum of the VOR signal shown in FIG. 3, a symmetrical image with respect to f.sub.e1 /2 and, after shifting to the left by 13280 Hz, the section of the spectrum to be processed, i.e. the frequency modulated carrier is centered about the frequency 3320 Hz. It is therefore noted that the new sampling frequency to be considered for the discriminator 5 is: f.sub.e =f.sub.e1 /2=13280 Hz, or four times 3320 Hz. In other words, the high pass output is in the required condition to be processed by the previously described discrimination procedure. More generally, using the processing system in FIG. 5, the reduced sampling frequency is deduced from the expression:

or:

The sampling frequency is therefore reduced to one third of its value with respect to the first solution, with everything happening as if the carrier frequency were itself divided by three. In this particular case, either the limited development of the expression (18) or the Arc cos function which slightly complicates the computations but provides greater accuracy may be used, equally, for the calculation of .OMEGA. and subsequently of y.sub.n.

The samples e.sub.n and subsequently x.sub.n (see FIG. 1) are then processed by a microprocessor, preferably a TMS 31010 made by the Texas Instruments Company, in order to carry out the computations mentioned above. The discriminator according to the invention is therefore entirely formed by a microprocessor, preferably provided with auxiliary memories, this microprocessor also being able to incorporate in the case of the VOR application at reduced sampling frequency, the band separator digital filter 7 and also, if appropriate, the digital filter 6.

The basic function of the digital discriminator 5 (see FIG. 1) is to determine the value y.sub.n by the expression (22) above, from a predetermined number N of consecutive samples x.sub.n taken between dates t.sub.q and t.sub.q+1, after which the next value of y.sub.n determined from the subsequent N samples x.sub.n is proceeded to (between dates t.sub.q+1 and t.sub.q+2). In order to avoid confusion in notation, it is recalled that strictly speaking, each sample x.sub.n corresponds to a distinct sample y.sub.n but by approximation it is considered, section by section, that the values of the N consecutive output samples y.sub.n are equal, this common value being precisely written y.sub.n for convenience of notation. N output samples can therefore be replaced each time by their representative y.sub.n which corresponds to an output sampling reduced by the ratio N with respect to the input giving an output frequency of f.sub.e /N. In the example described below the value of N is fixed at 64. Preferably, the sample x.sub.n are stored, in packets of 64, in an external RAM, hereafter referred to as the input-output port PA1 which forms means of accumulation of these input samples. Known buffer components and means of sequencing which are not shown enable the replacement of 64 input samples by the following 64 samples during a sequence of the computing program of y.sub.n following the reading of the samples, the processor working in real time.

FIG. 7 shows the general functional flowchart of the discriminator according to the invention; this flow chart, as well as the following ones, is formed from various boxes. In FIG. 7, box 101 indicates the start of the program. In box 102, the necessary initializations are carried out which consist, for the samples x.sub.n, in being set to read at the start of the memory area which contains them and, for the coefficients R.sub.01, R.sub.02, R.sub.11, R.sub.12 and R.sub.22, in loading their respective memory boxes referenced R01, R02, R11, R12 and R22 with the value 0. The first two samples x.sub.0 and x.sub.1 are read (boxes 103 and 104) and then, by means of a loop computation (boxes 105 and 106), the following sample is read and the coefficients R.sub.01 and R.sub.22 are computed. In order to do this, the auxiliary register ARO which was previously loaded with the value 63 is used as a counter, being decremented by one unit on each pass through the loop, until the value 0 is reached which marks the end of the computation of the coefficients and commands the exit from the loop. In box 107 it is determined which is the biggest of the two coefficients R.sub.11 and R.sub.22, after which (box 108) the other coefficients are put into the format of the value found in the previous box. The quantities R.sub.11 R.sub.22 -R.sub.12.sup.2, R.sub.02 R.sub.12 -R.sub.01 R.sub.22, R.sub.01 R.sub.12 -R.sub.02 R.sub.11 which respectively form the common denominator of the coefficients a.sub.1 and a.sub.2, the numerator of a.sub.1 and the numerator of a.sub.2 are then computed, in box 110 and compared with predetermined threshold values, which enables the values of a.sub.1 and a.sub.2 to be declared as valid or non-valid depending on whether the convergence criteria already mentioned above and materialized by these thresholds can be considered as satisfactory or not satisfactory. The computation of the current output sample y.sub.n (box 111) may then be proceeded to. Box 112 indicates the end of the general flowchart.

The flowchart in FIGS. 8a and 8b shows, between boxes 101 and 112, the more detailed progression of the flowchart in FIG. 7. Mathematical operations carried out in base two on the coefficients R.sub.11 and R.sub.22, as described below, enable these coefficients to be brought close to unity, without however standardizing them to 1, in order not to uselessly complicate the computations; consequently, it can be expected that the values determined subsequently for the coefficients a.sub.1 and a.sub.2 are not very close to 0 and 1 respectively. In practice, two positive threshold values are fixed s.sub.1 and s.sub.2 both close to 1, s.sub.1 by default and s.sub.2 by excess and it is required that, in order for the determination of the current output sample of the discriminator y.sub.n to be acceptable and declared valid, the values of the corresponding coefficients a.sub.1 and a.sub.2 must satisfy the following inequalities:

and

In the example described below, the values s.sub.1 =0.8 and s.sub.2 =1.2 are chosen. These two threshold values are in fact determined by experiment and can be changed at will in their respective 16-bit memory locations THRESHOLD 1 and THRESHOLD 2 during an operational calibration of the equipment. The decimal values 1.2 and 0.8 are produced from two respective values 1, contained in a memory location ONE and 3, which is 11 in binary, contained in a memory THREE. Firstly, the value 0.2 is determined by adding in the right hand section (fractional) of the accumulator the binary values 11, sucessively shifted towards the left by 4, 8 and 12 bits, which gives the binary value: 0.0011001100110011, a value very close to 0.2 in decimal, which is stored in the memory THRESHOLD 2 (box 121 in FIG. 8a). The value s.sub.2 =1.2 is in fact subsequently obtained by concatenation of the respective contents of the memories ONE and THRESHOLD 2 ((ONE)+(THRESHOLD 2)). Similarly, the value 0.8 is obtained (box 122) by adding the binary values 11 successively shifted by 2, 6, 10 and 14 bits and the result 0.1100110011001111 is put into the memory THRESHOLD 1. As mentioned above, the 64 samples x.sub.n are placed in order of arrival of an external RAM, at memory locations 147 to 210 of port PA1, which forms means of accumulation, the oldest sample being in 147. By means of another port PAO used as a counter incremented by unity on each IN instruction of the program, the samples x.sub.n are thus successively read and temporarily placed in a transit memory X2 of the internal RAM of the microprocessor. For more details on this particular possibility of memory extension, the Texas Instruments TMS 32010 User's guide 1983, pages 6-2 and 6-3, incorporated here for reference, can be referred to. In box 123, a program memory SERI is loaded with the value 147 and this value is then transferred to port PAO as a start value for incrementation. The memories allocated to the coefficients R.sub.01 to R.sub.22, namely R01L, R01H, R02L, R02H, R11L, R11H, R12L, R12H, R22L and R22H, which form the first means of storage, the letter L each time indicating the 16 bits of the low section and the letter H indicating the 16 bits of the high section of each coefficient are initialized, i.e. loaded with the value 0 (box 124). The register ARO used as a loop counter for the computation of coefficients R.sub.01 to R.sub.22 is then loaded with the value 63 (box 125), then the first two samples x.sub.0 and x.sub.1 are tranferred from the port PA1 to the memories X0 and X1 respectively (box 126). At the "START" label, the next sample, namely x.sub.2 is transferred into the memory X2 and the computation of the first term of each coefficient R.sub.01 to R.sub.22 according to the following expressions can then take place: ##STR1## In box 132 the following transfers between memories are carried out: (X1).fwdarw.X0 and (X.sub.2).fwdarw.X1. At 133, the auxiliary register ARO is decremented by unity and the new content of ARO (62) is tested in box 134. As this content is positive the "START" label is returned to, memory X2 is loaded with the next sample, namely x.sub.3 and the second term of the coefficients R.sub.01 to R.sub.22 is computed and added to the previous ones (boxes 127 to 131), according to the expressions 31 to 35. The computing loop is carried out as many times as necessary, the memories R01L, Ro1H to R22L, R22H contain the value of the coefficients R.sub.01 to R.sub.22 respectively, the register ARO then reaches the value 0, the test ARO=0 in box 134 proves positive and an exit from the loop takes place. For the computation of the coefficients R.sub.01 to R.sub.22, it will be noted that the samples x.sub.0 to x.sub.63 are, in their respective memories, written in 12 bits. Consequently, each term of x.sub.u x.sub.v type occupies a 24-bit format and the sum of 64 of these terms has a 30 bit format, which justifies the use of two memory locations (L and H) for the writing of each of the coefficients R.sub.01 to R.sub.22. In box 135 there is a test to find out which of the two coefficients R.sub.11 or R.sub.22 is the larger and, depending on whether it is R.sub.22 (box 136) or R.sub.11 (box 137), the program memory SERI is loaded with the high part of this largest value designated by SUP (R11, R22), i.e. with (R11H) or with (R22H). The continuation of the computation consists in loading the high part of the accumulator with (SERI) and in successively shifting this value by one bit to the left until a 1 appears in the MSB of the accumulator. If at the end of 15 successive shifts the MSB is still equal to 0, this indicates that (SERI)=0. It is then considered that the determination of y.sub.n is bad because it would be obtained from sample values x.sub.0 to x.sub. 63 which are too low, the computation of y.sub.n is not containued and the "ERROR" label is proceeded to directly (which is not shown on FIG. 8a). In the most probable case, the value 1 in the MSB is obtained at the end of a whole number r of shifts (r<15) and an exit from the computing loop takes place which operates the successive shifts at the rate of one per cycle. In order to carry out this computing loop (not shown) the decrementation of the auxiliary register AR1 may be used, for example, which is initially loaded with the value 15 and the auxiliary register ARO intended to register the value r after exiting from the loop. The procedure then goes on, in box 139, to another computing loop which serves to shift all the values of coefficients R.sub.01 to R.sub.22 by r bits to the left. To do this, the decrementation of the auxiliary register ARO is used which is initially loaded with the value r. At this stage of computation, the contents of the memories R01H, R02H, R12H and R22H which are considered to represent the whole values of the associated coefficients, are in a 16-bit format, the largest of them using the 16 bits. It is therefore possible to envisage completely ignoring the fractional parts of the coefficients R.sub.01 to R.sub.22 contained in the memories R01L to R22L respectively but, preferably, an additional rounding is carried out which consists in weighting by 1 positively and negatively respectively, depending on the positive or negative sign respectively of the corresponding coefficient, the most significant bit of the low part of this coefficient, after which only the new high part of the coefficient (after increasing or reducing by unity respectively) is retained for the rest of the calculations in order to simplify them (box 140). The computation continues (FIG. 8b) with the determination of the coefficients a.sub.1 and a.sub.2. Firstly, the common denominator of a.sub.1 and a.sub.2 is computed according to the expression: ##STR2##

Similarly, the respective numerators of the coefficients a.sub.1 and a.sub.2 which are stored in memories NA1 and NA2 respectively are then computed, according to the expressions: ##STR3##

For the three operations described above, it will be noted that the partial products of the type (R11H).times.(R22H) take, in register P and the accumulator of the microprocessor, a 32-bit format, as do the differences between partial products. In order to simplify the subsequent computations, only the content of the 16 bits of the high part of the accumulator (whole part) is stored in the memory after possible rounding, for the common denominator, while the respective values of the numerators of a.sub.1 and a.sub.2 are each stored in two memory locations (whole part and fractional part). The value of the coefficient a.sub.2 is computed by division: ##EQU27## The result of the division which uses 32-bit dividends is stored in two memories associated with the division subroutine DIV32, QUOTIM (whole part) and QUOTIL (frational part). After this it is first checked that the result is positive using the test (QUOTIM)>0, which is implicit by virtue of the above inequality (26). If this is not the case the "ERROR" label is used. For (QUOTIM)>0, the convergence test of the inequality (26) is carried out (box 145) i.e. the following operations are carried out successively:

If the result of this operation is negative, the "ERROR" label is used, otherwise the test is continued with:

If the result of this operation is positive, the "ERROR" label is used, otherwise the computation of a.sub.1 is effected: ##EQU28##

The result of the division is stored as before in the QUOTIM and QUOTIL memories (in the place of a.sub.2) which form second means of storage. In order for the value of a.sub.1 to be acceptable, i.e. sufficiently close to 0 by positive or negative values, it is checked that it satisfies the inequality test (27) (box 147), i.e. the sign of a.sub.1 is first tested:

if (QUOTIM)+(QUOTIL)>0 the operation is carried out:

if this result is positive "ERROR" is used, otherwise ((QUOTIM)+(QUOTIL)>0), the operation is carried out:

if this result is negative, "ERROR" is used, otherwise the computation of a.sub.1 is declared valid and the computation of y.sub.n can then be effected.

In the present example of embodiment it is preferred to compute the value of y.sub.n by application of the above formula (22)', giving, using the notations chosen for the description of the flowchart in FIG. 8: ##EQU29## where ##EQU30## K.sub.1 is a multiplicative constant that can generally be omitted. In fact, the ultimate aim of the discrimination carried out on the composite VOR signal which constitutes the preferred application of the invention consists in the real time determination of the phase called the reference phase of this signal. This purpose is achieved by isolating a 30 Hz frequency line, this single sinusoid being obtained in the form of y.sub.n ordinate points provided at regular intervals. From the values of y.sub.n it is then possible, using a digital method, an example of which is given below, to extract the amplitude and phase of the sought line. Now, the coefficient K.sub.1 only affects the value of the amplitude while it is basically the phase which is sought. Also, it will be noted that the correct value of the amplitude can be obtained at the end of the computations by bringing the multiplicative coefficient K.sub.1 into effect only at this stage. The value y'.sub.n representing y.sub.n actually calculated is therefore: ##EQU31##

For this, the following operation is first carried out

Given the format (NA2) which is in fact a 32-bit format, in the memories NA2H and NA2L, rounding is preferably first carried out, by adding 1 just after the point, after which only the 16 MSB are stored in memory NA2H and it is this latter value which is used for the multiplication mentioned above. The result is supplied in 32 bits to the register P. A subroutine RACAR is then called which carries out the operation:

X being a memory location which therefore contains 16 bits.

After this, NA1, is transferred from memory locations DLOWM and DLOWL to the accumulator and the final operation indicated by expression (44) is carried out by calling a subroutine DIV 32. The value of y'.sub.n thus obtained is stored in the form of a whole part in 16 bits and a fractional part in 16 bits in memories QUOTIM and QUOTIL respectively which form the second means of storage.

The digital method used to extract the phase of the signal at 30 Hz is, preferably, the method called the "recursive discrete Fourier transform" or "recursive DFT" method which is explained in its general principles in the appendix. The problem raised reduces to computing: ##EQU32## Y(k) being a complex number which represents in amplitude (a) and in phase (.phi.) the 30 Hz spectral line which is in question.

In the expression (45) :

M is the whole number of distinct successive values y.sub.n chosen for the computation, representing several sinusoids of the modulating 30 Hz signal, for example M=83 (12 sinusoids for a frequency f.sub.e /N of 207.5 Hz of output samples (y.sub.n);

k is such that ##EQU33## giving, with the previously chosen values: ##EQU34## k=12

It is desired to determine Arg Y(k) and .vertline.Y(k).vertline.. To do this use may be made of the fast algorithm given in the appendix and defined by:

v.sub.-1 =0

v.sub.-2 =0

for

q=0 to M-1 ##EQU35## after M values (real) of v.sub.q have been computed: ##EQU36## For each sample y.sub.n (i.e. y.sub.q when working with respect to M) at the output of the discriminator 5, the computation of the corresponding v.sub.q according to expression (47), in the considered section 0 to M-1, therfore takes place:

In practice, in order to be able to double the refresh rate of the memories, it is preferred to interlace by 50% the determination of Y(k), i.e. for each values of y.sub.q two values of v.sub.q are computed, a first v.sub.q value resulting from a pair of earlier values V.sub.q-1, v.sub.q-2, previously calculated and respectively contained in the memories V1 and V2 and a second v.sub.q value resulting from another pair of earlier different values v.sub.q-1, v.sub.q-2 contained in memories V10 and V20 respectively. The interlacing is such that, when the position is maximum for a pair of values v.sub.q-1, v.sub.q-2 or M-1 and M-2 (82 and 81) and is suitable for the associated value of Y(k) to be computed according to formula (48), the position of the other pair of values v.sub.q-1, v.sub.q-2 is equal to M-1/2 and M-3/2 (41 and 40) for the other pair of values v.sub.q-1 .multidot.v.sub.q-2.

The ports PA1 and PA0, respectively working as an external RAM and as a counter with initial setting and associated incrementation, can again be used for placing the values v.sub.q-1, v.sub.q-2 and v.sub.q necessary for the computations into memory.

The flowchart in FIG. 9 shows the computation of the two v.sub.q values associated with each y.sub.q value. Box 151 indicates the start of the program which, preferably, follows the program of computing y.sub.n (FIGS. 8a and 8b). In box 152 the value of 2 cos 2k.pi./M is read. This value, which for example is stored in memory location 290 of PA1 is transferred, by means of PA0 to location RXN9H of the internal data memory of the TMS 32010 microprocessor. In box 153 the contents of memory locations 608, 609, 610, 611 are transferred in a similar way by a series of IN instructions into memories V1, V10, V2, V20 respectively of the microprocessor. Box 154 symbolizes the following operations for the so-called recent DFT computing sequence which concerns (V1) and (V2):

(V1).times.(RXN9H).fwdarw.ACC

(ACC)+(QUOTIM).fwdarw.ACC

(ACC)-(V2).fwdarw.ACC

(ACC).fwdarw.V2

which carries out the literal operation:

or the operation: ##EQU37## The new content of V2 thus changes from the value v.sub.q-2 to the value v.sub.q (V1) remaining unchanged at this stage. Box 155 symbolizes the same operation as in box 154 but for the so-called more than 50% DFT computation sequence which concerns V10 and V20. The new content of V20 thus changes from the value v.sub.q-2 to the value v.sub.q, (V10) remaining unchanged at this stage.

In box 156, PA0 is loaded with the value 608 and then, by a series of OUT instructions, the contents of memory locations V2, V20, V1 and V10 respectively are transferred in this order into memory locations 608, 609, 610 and 611 respectively of port PA1.

By means of this latter storage operation the necessary substitutions are made in order to change from the values v.sub.q, v.sub.q-1 to the values v.sub.q-1, v.sub.q-2, for each pair of adjacent values v, which is necessary for the purpose of the following computation (FIGS. 8a, 8b and 9) which is precisely characterized by an increase by 1 of the mute indices n (FIGS. 8a, 8b) and q (FIG. 9). Box 157 marks the end of the computation and the return to box 123 (FIG. 8a) for the computation of the following values of y.sub.n (y.sub.q), v.sub.q and v.sub.q-1. A count, which is not explained here, indicates which pairs v.sub.q-1, v.sub.q-2, are identified with the pairs v.sub.M-1 and v.sub.M-2 from which the values of Y(k) are computed each time, according to expression (48), Y(k) being a complex number from which can be extracted the modulus and particularly the argument .phi. which is extremely helpful to known. This computation, which is specialist in the field is capable of programming, is carried out by the TMS 32010 microprocessor and each value of .phi., called .phi..sub.ref is then transferred by the data bus according to a certain protocol to a FIFO memory forming the second means of storage, called the results FIFO memory. By connecting, (in a way that is not shown) an instruction FIFO memory arranged in parallel with respect to the results FIFO memory on the data bus, a buffer memory element can be created between the TMS 32010 used as a slave microprocessor and another microprocessor used as a master microprocessor, a MOTOROLA 6809 for example, the latter synchronized with the TMS 32010, capable of using the results contained in the results FIFO memory and also carrying out management tasks which are outside the scope of the invention. The information available at the "ERROR" label is transmitted by the results FIFO to the master microprocessor which then decides upon the possible modification of the values of the s.sub.1 or s.sub.2 thresholds or the increased amplification to be given to the analog input signal or a predetermined delay to be observed before repeating the computations which the discriminator 5 must carry out.

It has been seen above, with reference to FIG. 1, that in the presence of an input signal having too much noise it was advantageous to introduce a pre-filtering by means of an input filter 6. The inclusion of the filter 6 in the circuit can be controlled by an appropriate signal on the conductor 4 from the TMS 32010 microprocessor itself, which is represented in FIG. 1 by the broken line 11; instead of the value of the S/N ratio a value representing the noise power .sigma..sup.2 can be used and the switch 3 can be activated in the position shown in FIG. 1 when this representative value exceeds a certain threshold S" to be determined during the equipment test. As a representative value of .sigma..sup.2, expression (24) is chosen:

The coefficient R.sub.11 or the coefficient R.sub.22 can be equally used for R.sub.0 and, in order to carry out the computation according to expression (49), it is previously necessary to store the quantities a.sub.1 and a.sub.2 in their own memory locations in the microprocessor.

For the final computation of the phase of the signal to be discriminated, numerical methods other than the recursive DFT can be used, such as a method of least squares for example.

In its most general form, expression (1) is written:

where

t: time of arrival of the sample

c: d.c. component

b.sub.t : zero average additive noise

In vector form, the expression (50) is written ##EQU38## or again:

where:

H.sub.t =[1 sin .omega.t cos .omega.t] and

X=[c a cos .phi.a sin .phi.].sup.T

The problem consists in estimating the quantities c, a and .phi., i.e. the vector X assuming stationary or pseudo-stationary conditions.

Multiplying both sides of equation (51) by H.sub.t.sup.T there is obtained:

Accomulating over M points: ##EQU39## the noise b.sub.n being of zero average and not correlated with the terms in sin t and cos t, can be ignored and the result is: ##EQU40## being a symmetrical square matrix, and:

where ##EQU41##

a and .phi. can then be extracted by calculating P and Z and by solving the matrix equation:

It will be noted that, assuming an equidistribution or regular distribution of observation times t.sub.q and observation duration T.sub.1 that is a multiple of the analysed signal period T(T=2.pi./.omega.), it can be shown that the matrix P tends towards a diagonal matrix: ##EQU42##

Assuming the interval between samples as being constant, for M regularly spaced points during the interval T.sub.1, the expression of sampling times t.sub.1 becomes:

The Z column vector is then written: ##EQU43## As the duration of observation is a multiple of the period of the component which it is desired to estimate,

can be written,

k being an integral multiplicative factor

Z is then written: ##EQU44## PS Y(k) designates the DFT (Discrete Fourier Transform) of the series y((q T.sub.1 /M) for values of q from 0 to M-1.

Expression (52) becomes: ##EQU45##

The computation of X therefore can be summed up in the knowledge of the lines of order 0 and k in the DFT of the series y(q T.sub.1 /M. In order to find the complex number Y(k) it is preferable to use a recursive DFT method. It will be noted that knowledge of Y(k) is sufficient for estimating the vector X. In fact,

and

can be written, which means to the nearest multiplicative factor 2/M there is:

and

The problem can be summed up by estimating Y(k) such that: ##EQU46## To facilitate the writing, this is put as:

It is noted that Y(k) is a polynomial in powers of W whose coefficients are the samples of the received signal.

Y(k) is broken down in order to make the various samples in the series y(0), y(1), . . . , y(q), . . . , y(M-1) appear. This is put as:

(from which: W=e.sup.-j.theta.

With all computations complete, the following is deduced from them:

where:

v.sub.-1 =0 and v.sub.-2 =0

Y(k)=e.sup.-j.theta. v.sub.M-1 -v.sub.M-2

Finally the following algorithm is obtained: initialization:

on each received sample y(q):

v.sub.q =y(q)+2 cos (2k.pi./M)v.sub.q-1 -v.sub.q-2 (real value)

at the end of the sequence of observations:

Y(k)=e.sup.+j2 .pi./M v.sub.M-1 -v.sub.M-1 (complex value)

It will be noted that the reduction in the number of operations with respect to the standard DFT strengthens the possibility of a so-called on-line computation, i.e. inserted in the time intervals between samples. This enables the value of y(k) associated with the M previous samples to be obtained shortly after reception of the last sample (or order M-1). Also, this on-line computation is accompanied by a significant reduction in the memories necessary for the computation. This reduction is even greater in this case than in other recursive methods because of the small number of variables.