1. Preliminary Overview of the Raw Data
In the westside of the La Spezia port—ITA (φ = 44 03 59.58N, λ = 09 50 49.22E, elevation e = −4.5 [m], date 14.09.2012, time Δt 09.10am: 09.20am (GMT)) the magnetogram of Figure 1 has been observed by means a magnetovariograph (sampling rate s = 5 × 10−1 [sec]). The time variation of the geomagnetic field vertical component ΔZ is studied in the present work. The magnetovariogram (Figure 1) consists of 1200 samples for a length L = 10 [min]; ΔZ has an accuracy of 1 × 10−1 [nT] and is returned in the graphs (and in the calculation) with precision of 1 [nT] (approximation by truncation). In theory, the measured function solves a spectral window Λmin = 2 [s] → Λmax = 300 [s]. The planetary index K of geomagnetic field activity during the experiment was = 1 (standard K → H comp.).
Preliminary observation of raw data shows a typical harbour high-noise coastal proximity field. It is characterized by a very large spectral window dominated by a high frequency noise band close to the Nyquist frequency (which probably carries important aliasing components foreign to the object of the present study) interfered by a lower frequency noise band. The magnetovariogram qualitative observation (observed field) highlights some preliminary physical informations on the origin of the local electromagnetic noise. For description of the parameters of our interest, we divide this magnetogram in 4 windows (W1, 2, 3, 4) plus a continuous subset (W5) of the window 4. The magnetogram is characterized by a very large spectral window dominated by a high frequency noise band close to the Nyquist frequency (which probably carries important aliasing components foreign to the object of the present study) interfered by a lower frequency noise band (Figure 1).
- W1 (samples 1 → 102): particularly disturbed field due to the unfolding of the chain of magnetometers. One of these devices provided the magnetogram di.
- W2 (samples 122 → 149): noise generated by the power supply test. In correspondence with the shift AC → DC there is a intensity drop of the ΔZ signal background of about 20 [nT]. The phenomenon is related to the electromagnetic activity of the generator and its geometric position with respect to the sensor considered.
Figure 1. Underwater magnetic field in harbour of La Spezia (ITA). ΔZ instrumental sensibility 0.1 [nT]. Wn: subsets of qualitative analysis of the magnetic field (Raw data). Unit of measure of the magnetograms: X = [sec • 0.5], Y = 1 [nT]. (These dimensions are valid for all magnetograms of the paper).
- W3 (samples 150 → 298): transition to AC power supply with increased intensity of the induced magnetic background (ΔZ) of about 20 [nT]. Time variation dZt/dt
Probably dZt/dt is not related of the magnetic noise produced by the electrical supply system but it is coming from an unknown artificial source .
- W4 (samples 298 → 1201): definitive return to the power supply by generator and inversion of the phenomenon of magnetic interference (see W2) on the background of the ΔZ signal.
- W5 (samples 750 → 1050, sub-window W4); decrease of the high frequency noise (noise W5 ≃ 20% noise W4). The phenomenon is probably dues to the temporary shutdown of a high frequency noise source close to the sensors. In the qualitative discussion of the characteristics of the total magnetogram, this fact requires that the subset W5 is defined individually and not as a part of W4 that contains it. On the other hand, the previous and next data set W4 is more homogeneous and for this reason it is preferable to define W5 as a subset of W4 rather than dividing the data into three subsets.
Furthermore, the magnetogram is characterized by the presence of probably instrumental spikes (same times classifieds as “electronic disturbance”). We observe also several electronic spikes and an artificial edge imposing two numerical interventions to stabilize the series (cleaning data action): 1—to delete spikes and 2—to delete edge induced by the AC/DC power passage and vice-versa. The standard information processing protocols on data whit high time transient would require the application of numerical techniques for signal strengthening (i.e. FFT filtering) to make the signal information more usable. In the present case we do not apply these procedures to subject the synthetic signal reconstruction protocol to maximum operative stress.
2. Data Clearance
The data cleaning action is developed in two distinct steps: first one elimination of electronic noise (spikes), second one elimination of the edge (man made noise, data 291 - 292).
2.1. Killing the Spikes
We use a standard numerical technique  . It is is a very wide spectrum low pass filter that excludes impulsive signals of very high frequency. At this action both the spectrum of the CMT and that of the magnetic activity of the target (divers) survives. In the present example are classified and erased the electronic spikes. These signals have wavelength Ʌ ≤ 2Ñmeasure = 1.0 [sec] (Ñmeasure = sampling rate = 0.5 [sec]). The cutoff wavelength applied for our LP filter is Ʌcut = 1.5 [sec] after an oversampling operation performed to stabilize the action of the Fourier procedure for the frequency band of our interest . In Figure 2(a) is shown the result of the LP filter action: the magnetogram has not electronic spikes.
Figure 2. Underwater magnetic field in harbour of La Spezia (ITA). ΔZ field cleaned of electronic spikes (a). Amplitude step in the magnetic field Z component; it is originated by human activity (b).
2.2. Killing the Sharp Edge Step (A291 → 292)
The series of raw data there is an artificial singularity generated by the change in the power supply of the measurement network from which the data of the paper is coming from. Between samples 291 and 292 there is a signal drop of about 20 [nT] (Figure 2(b)). The Fourier procedure has a great difficulty in describing this type of signal geometry which we can define as a cutting pitch (transition 291 - 292). This geometry generates pollution problems of the information disseminated to the entire series of recorded data (particularly annoying in the band of high frequency).
These problems can be referred also to the well-known phenomenon of aliasing. To remove this problem (after elimination of the spikes) we extract a subset of 110 samples centered on step 291 - 292 (Figure 3).
In this series the discontinuity is between the sample 51 and the sample 52. To obtain a value of discontinuity representative of the series and not only of the two samples 51 - 52, we assign to the sample 51 the value of the linear regression of the interval 1 - 51 (series 1) and to the sample 52 that of field 52 - 110 (series 2). The difference between these two values is adopted as compensation “k” factor (Figure 4). We define K as series equalizer factor.
The equalization of the total series is obtained by the relationship Series (1 - 51) + k = Series (1 - 51) + (−20 [nT]) (Figure 5). The “equalized” series is shown in Figure 6, it is base data set of our job. We define this data series as “observed field” function (Fobs).
Figure 3. Elimination of electronic spikes and extraction of subsets (series 1 and 2) for the contrast of the singularity (man made step): the general flow of calculation.
Figure 4. Computation of the K “equalization factor”: the general flow of calculation.
Figure 5. Elimination of the singularity (man made step) by means K “equalization factor”: the general flow of calculation.
Figure 6. Cleaned and equalized magnetograms. Observed field function Fobs ΔZ.
3. The Signal Synthetic REConstruction (SYRec): A New Choice for the Low Interference Signal Reconstruction
As is well known, the discontinuities of data in the numerical series generates falls (loses) of direct and indirect information . We define information directly lost the ones contained in the time of unregistration and indirectly lost (sometimes called the second order information fall) those are due to the spectral pollution coming from the interruption. Longer is the crash window (interruption of the data sequence) higher is this pollution level. The main cause of this action is the practical impossibility of the Fourier calculation protocols to define the geometric conditions of non-derivability (sharp edges).
In general the aim of the reconstruction of missing data is not to retrieve the information carried by the lost data but to protect the information carried by the recorded data. Numerical reconstruction techniques are many, their spectral effects are often related to the ratio transience of the sampled function-sampling density . The most satisfactory reconstructive solutions are heavy, require not short calculation time and therefore are not compatible with control systems that base their effectiveness on the immediacy of the response . Moreover, if the systems are made up of many peripherals (e.g. a magnetometer array that observes the magnetic singularity variations over a large area), it is essential to delegate part of the measurement intelligence to each peripheral so as not to compromise the effectiveness of the system in any eventual single instrumental crash For this reason, heavy computational protocols are not operative in the case, for example, of environmental magnetic control-reaction systems. But on the other hand it is fundamental to equip the system with a numerical suture protocol for unregistered periods. In fact, the crash condition of duration T can generate a much longer system blindness condition. In the present case a fall of the system of duration T = 120 [s] samples = 60 [sec] blinds the system for a period T = 3 [min]. More, operational practice shows that pollution diffusion is greater and more extensive on harmonics of our interest (higher frequency harmonics) . If the system must act on these harmonics it is necessary to proceed to a suture as effective and fast as possible. Working towards the high frequencies (relative to this type of signals essentially expressed in the low frequency band) we have the following problem: an effective suture is necessary to contain the pollution of the information and at the same time it is necessary to use light and fast calculations for not compromise the response speed of the system. The solution adopted in the magnetic control system of the HaPS Harbor Protection System (EDA) project is the “SYRec” algorithm (SYNTHETIC REConstruction). SYRec calculates a sub-sampled function (LFf Low Frequency function) with sampling rate = ÑDIRAC = Δtcrash = Lcrash and phase equal Δtcrash. From LFf we extract a series of lengths Δtcrash and phase Φ = ΦΔtcrash and then insert it in the observed series crash window. This step is the reconstruction of the low frequency component. To this component is added (in inverse progression) the signal of length Δtcrash measured immediately before the crash (HFf High Frequency function). The response of the merge is a numerical series of length L = Δtcrash in phase with Δtcrash containing the low frequencies of the entire recording and the high frequencies closer to the crash period.
We discuss and compare the results obtained by the SYRec standard with those of two standard fast suture actions:
- SYRec standard (signal SYnthetic REConstruction). Enough fast and highly effective in containing information pollution.
- LFR low frequency reconstruction. Enough fast but unsatisfactory for high frequencies.
- HFR high frequency reconstruction. Fast but unsatisfactory (in some spectral conditions harmful).
To start the study of the effectiveness of the suture protocols we generate the interrupted magnetogram of Figure 7 from the observed magnetogram of Figure 6 (noise condition = hard) and proceed to the synthetic reconstructions. The information contents of the reconstructed magnetograms coming from LFR, HFR, SYRec techniques will be compared with those of the continuous (and clean) observed magnetogram (Figure 6).
We start generating the broken magnetogram with an artificial interruption of length Δt = 30 [sec] between the sample 452 and the sample 512 (border amplitudes values of the crash period = 48 [nT], 54 [nT]) (Figure 7). The artificial interruption is long ΔT
3.1. Low Frequency Reconstruction—LFR
The LFR procedure starts with the extraction of a continuous series of data from
Figure 7. Artificial interruption of length Δt = 30 [sec], sample 452 → 512. BcW and AcW: computation data windows used to the the reconstruction of the interruption (data series “A” and “B”). Above the “Dirac Comb” for the subsampling for the construction of the low frequency function.
the interrupted series (observed Fobs function) for a subsampling action execute with a sampling rate equal to the length of the signal interruption
where signal interruption end time
and signal interruption start time
and subsampling action b in phase with the interruption (Figure 8).
This action produces a continuous subsampled function (Fus). Then we proceed to a resampling of Fus by polynomial approximation of 5 order with a sampling rate Ñpol equal to the Fobs one.
The result function (called the Fpol polynomial function) contains the low-frequency information of Fobs (except for computational approximations) and it has sampling rate and length equal to Fobs.
le due funzioni hanno uguale densità numerica ρ e fase φ e sono quindi confrontabili in TD (time domain)
Obviously the lower the length of the signal interruption the better the reconstruction in accordance with the general rule
The resampling action has a cost in terms of computation weight. This cost is directly proportional to the length of the sector subjected to resampling for this reason it is necessary to make a compromise in the choice of this length. In the present case we propose a calculation segment of length Lpol
where BcW (Before crash Window)
and AcW (After crash Window)
This choice, in our opinion, is the best compromise between the weight of data processing and effectiveness in limiting the numerical pollution due to the suture of the crash.
Finally, the 451-→ 513 data window is extracted from the Fpol and it is inserted in phase in the crash window of the Fobs series.
Where the sample 451 is time of the start of the crash ti and 513 is the time of the end of the crash tf.
The graphic performance of the LFR reconstruction is shown in Figure 9.
This signal reconstruction technique is excellent for protecting the physical information of medium-low frequency signals less for those of high frequency
Figure 8. Low frequency information extraction. Ñ of the Dirac Comb = Δt crash (interruption). The sampling rate is equal at the interruption and in phase with this one. The sub-sampled function is continuous.
Figure 9. Low frequency reconstructed “LFR” signal.
signals (particularly if Lcrash is long). It is not a heavy numerical technique but it is too penalizing for high frequency observations.
3.2. High Frequency Reconstruction—HFR
HFR is a very simple and fast calculation option. It is based on the action of a counter that detects on quasi-real time the absence of the data in the measured serie and replaces them with the corresponding series immediately preceding the metrological crash (in reverse sequence). For example if the values n + 1, n + 2, n + 3, n + 4 are lost they are substituted with the data sequence is n, n − 1, n − 2, n − 3 (Figure 10).
The operation stops when the counter detects a new measured data. In a qualitative way we can affirm HFR technique is reliable if the series in question is stationary (or with a low temporal increase) and if its spectral window is not too large . Otherwise HFR is not reliable as it can create pollution on the spectral informations more or less as much as the crash effects and even more in borderline cases. We proceed to the reconstruct of the interruption Δtcrash = t/10 by means the HFR technique. Figure 11 shows the effect of the present HFR application (not-stationary series).
Figure 10. Action of the data acquisition counter for the reconstruction of the HFR high frequency signal.
Figure 11. High frequency reconstructed “HFR” signal.
Figure 12(a) shows the spectral effects of HFR in the case of a stationary signal, in particular the black line is the amplitude spectrum of the original series (continuous), red line is the reconstructed HFR one. In Figure 12(b) the HFR is applied to a non-stationary series. The inconsistency of the HFR spectrum (red line) with respect to real one (black line) is clearly observable also in qualitative way.
The percentage error “e%” of HFR is defined in the graph of Figure 13. In this figure the black line represents the error in case of stationarity while the red one in case of non-stationarity.
The percentage error produced by HFR in a stationary numerical environment can be considered acceptable while that for non-stationary numerical environment is not acceptable.
If the original numerical series is not stationary, the HFR has no acceptable capacity to protect spectral information.
3.3. Synthetic Reconstruction—SYRec
As seen reconstructions of the LFR and HFR signal show advantages in execution speed but also heavy disadvantages in information protection performance. In particular, LFR has a calculation speed compatible with an environmental
Figure 12. Comparison between the continuous signal spectrum (Fobs) and the HFR signal spectrum. Black line Fobs, red line HFR. Graph “a” HFR reconstruction of an almost stationary function; graph “b” reconstruction of a non-stationary function. X order number, Y spectral amplitude.
Figure 13. Comparison of the informative efficacy of HFR in stationary and non-stationary function conditions. Black pointed line difference between the HRF spectrum and continuous stationary function spectrum, Red dotted and pointed line in non-stationary continuous function. X order number, Y difference in percentage energy “e%” for some elementary signal.
control-reaction system but does not sufficiently protect the high frequency band (our maximum interest). While HFR is very fast but very dangerous in the case of non-stationarity of the measured function (Fobs). Reconstruction (SYRec) procedure is the merge of the numerical actions of LFR and HFR. It is built to merge the best qualities of the two standard techniques of reconstructions. SYRec sums in Time Domain (in phase) the series of data of length Δtcrash coming from LFR and the correspondents ones of HFR. In this way SYRec integrates the control of low frequencies (from LRF capability) to that of high frequencies (from HFR capability).
According to Figure 14 we propose the following calculation flow consisting of three fundamentals steps:
Figure 14. Calculation flow of the synthetic reconstruction algorithm “SYRec” of the high noise and not continued signal.
- Definition of the number and position of lost samples (the subset of Fobs of length Lcrash = tcrash) by means a sequential counter.
- Lost data reconstruction by means the merge of the low and high frequency information according to LFR and HFR procedures.
- and then to insert on the in the Lcrash window of Fobs this composed data series.
The sub-sampled series Fus is extracted (in phase) from Fobs crashed by means of a Dirac subsampling function (Figure 8) where:
The undersampled series (Fus) is continuous and transparent to the interruption Δtcrash window and it has
From Fus the polynomial function Fpol of length 4Lcrash is extracted (in the present case we decide 4Lcrash to be best length for information protection effectiveness ratio/calculation weight) and sample step Δtpol = Δtobs
and sampling rate
from the Fpol we extract the LFR subset data of length Lcrash in phase with Δtcrash (Figure 8) with the following reading function WLFR
where ti is the crash window start time and
where tf is the crash window end time (Figure 15).
The high frequency component is added to the low frequency over-sampled serie with reverse phase. The final result of the action is shown in Figure 16. The computation time paid for this action is not too high (HFR computation time + LFR computation time = SYRec computation time). HFR computation time + LFR computation time is not much higher than the LFR calculation time. In general, the slowdown in SYRec production is less than 10% of the LFR calculation time.
4. Spectral Comparison of the Reconstruction Algorithms Effectiveness
To quantify the efficacy of LFR, HFR and SYRec in information protection we compare the amplitude spectrum of the original causal function Fobs cleaned and continuous with the spectra of the three continue reconstructed causal functions F(LFR), F(HFR), F(SYRec) Figure 17.
To protect this action from edge effects all the series considered are previously subjected to smoothing by means the smoothing function so called “cosine bell” (Figure 18).
The action of the cosine bell function is defined in (18) where Fn is the result of the cosine bell smoothing
This datum manipulation produces spectral stability benefits unrelated to the type of suture technique adopted. This increase of spectral stability is about the same for the three techniques of reconstruction (LFR, HFR, SYRec) and therefore their spectral comparative analysis does not lose validity.
We observe that in the frequency domain the cosine bell produces heavier distortive effects especially when applied to short numerical series.
The quantification of these distortions is easily computable from (19) (FFT of cosine bell)  :
Figure 15. Numerical parameters of the box car subject of the reconstruction in the present paper. X = amplitude, Y = number of samples. L box car length.
Figure 16. Synthetic reconstructed “SYRec” signal.
Figure 17. Amplitude spectra of the continuous signal Fobs and of the continuous signals reconstructed by means of the HFR, LFR and SYRec protocols. X axis Amplitude, Y axis number of order.
Figure 18. Data manipulation: box car parameters used to counteract the border differential effect in the spectral stability comparison between the continuous signal and the reconstructed signals L = 960 samples, M = 120 samples, N = sample, Ñ = 0.5 [s], A = A[nT].
but also these effects are distributed in approximately the same way on the Fourier Transforms of the Fobs, FLFR, FHFR, FSYRec and therefore do not intervene in the effectiveness of their comparison.
This condition justifies the validity of the differential spectral comparison of WLFR, WHFR, WSYRec with Wobs where
Figure 19 shows the spectral difference between Wobs and repetitively WLFR, WHFR, WSYRec.
The observation of the spectral difference indicates:
- the reconstruction of the signal performed in HFR pollutes both the high frequency band and the low band and therefore it is not suitable for an accurate reconstruction respecting the information capacity in the signal itself;
- the LFR reconstruction solves, in large part, the HFR problems in low frequency band but is more or less transparent to high frequency band and therefore its performance is not sufficient (especially for high frequency studies);
- SYRec contrasts in a very effective manner the spectral distortions both for the low and high frequency components. The price paid is increased computation time. But this cost does not exceed (in general) about 10% of the computational time of the other two data reconstruction procedures. It is acceptable.
Figure 19. Spectral difference between the Fobs continuous function and the corresponding reconstituted HFR, LFR and SYRec functions.
5. The SYRec Standard in Observatories Activity
To complete the comparative analysis of SYRec spectral efficacy, we compare its results with those of the low-frequency reconstructed “borders smoothed” LFR_bs (Figure 20).
This suture technique is able to contain spectral instabilities of both high frequency and low frequency but is very heavy as calculation time. It is therefore valid for “data collection” and observatory studies  but it is not suitable for “quick reading-response” activities of environmental control systems. This comparison allows us to test the efficacy of SYRec with a synthetic high stability signal reconstruction procedure.
According to Figure 22 in order to obtain the LFRs suture we proceed to the calculation sub-sets A and B
The results of the smoothing process are shown in Figure 23.
We compare at last the efficacy of SYRec action with LFRs suture protocol in use in our reference magnetometric stations (Figure 24). The gain in spectral stability obtained by LFRs (Figure 24(e_LFRs)) is qualitatively compared to that of SYRec (Figure 24(SYREC)).
The difference in deviation between the two reconstruction methods spectra with respect to the original continuous function spectrum is evidently greater for
Figure 20. Subsets “a” and “b” to be subjected to the “border smooth BS” computation procedure to obtain a LFRs type suture.
Figure 21. Graphic expansion of the subsets “a” and “b” submitted to border smooth computation (from Figure 20).
Figure 22. Data manipulation: quantitative parameters of the “border smooth” protocol action, X = sample, Y = amplitude [nT].
Figure 23. Border smoothing effect in subsets “A” and “B” of the LFR_BS algorithm.
the LFRs technique. This fact is fundamental for the study of the real energy distribution between the elementary harmonics of the reconstituted function and therefore fundamental for magnetic measurements of singularity  . SYRec is preferable to LFRs both in terms of calculation speed and performance stability of results. Following to this performance, the SYRec protocol has been included (dec. 01 2018) also in medium-term magnetic observation stations
Figure 24. Spectral difference between the Fobs continuous function and the corresponding reconstituted LFRs.
(MTM) (4 - 12 [days]) (magnetoimetric reference stations for detection of low energy magnetic signals, from quasi-point sources and kinetics active in high electromagnetic noise environment).
The purpose of SYRec protocol is the production of a suture system of interrupted magnetic recordings not very demanding for the automatic calculation area and highly effective in protecting the physical information of the recording. This target is vital in operative evaluation of detection systems of low energy, quasi-punctiform and kinetic sources in high noise magnetic environment. The performance of the SYRec numerical protocol was tested in over 80 trial actions performed in critical high noise conditions both in an underwater environment (protection of ports) and in terrestrial environments (protection of critical and urban structures), both military and civil. The SYRec response was compared with the most commonly used fast techniques to suture of interrupted numerical series (LFR, HFR). This operational has shown a decisive improvement in the protection of information. This result is obtained by paying an increase of the calculation weight not exceeding 10% respect to LFR techniques (more effective than the HDR technique but also more demanding from the computational point of view). This information security protocol worked for over 240 consecutive hours (Isola Rossa observatory) without showing instability. In this context, a very clear gain in terms of computation time and weight was also highlighted. The SYRec computational procedure was therefore included in the medium-term-magnetometric observatories (MTM). Today, this procedure validated and its difficulty were structuring started.
This study is part of EU and Italian Harbor/Base Protection research action. It was supported by SEGREDIFESA of the Italian Ministry of Defence under the National Military Research Plan R & T, projects C.A.I.Ma.N., La.Ma.1.0_2.0 and European Defence Agency by Project Ha.P.S. (SWE Lead Nation, ITA, GER, NOR). Field trials were managed by MARI.CO.DRAG Marina Militare in La Spezia Harbour (ITA) and by WTD71 Bundesmarine in Eckerfoerde Horbour (GER). This study was coordinated by CSSN-MM. Isola Rossa Geomagnetic Base Station (ITA) was managed in logistic collaboration with the S&T CMRE SP (ITA) (Capo Teulada Experiment-2007).
Thanks to all of you.
 Faggioni, O., Soldani, M., Gabellone, G., Hollett, R.D. and Kessel, R.T. (2010) Undersea Harbor Defence: A New Choice in Magnetic Networks. Journal of Applied Geophysics, 72, 46-56.
 Wax, M. and Kailath, T. (1985) Detection of Signals by Information Theoretic Criteria. IEEE Transactions on Acoustics, Speech, and Signal Processing, 33, 387-392.
 Faggioni, O., Soldani, M., Cozzani, G. and Zunino, R. (2018) Informative Signal Analysis: Metrology of the Geomagnetic Singularities in Low-Density Ionic Solution (Seawater). Journal of Signal and Information Processing, 9, 1-23.
 Faggioni, O., Palangio, P. and Pinna, E. (1991) Geomagnetic Observatory “Terranova Bay Station-Antarctica”: Synthetic Reconstruction of the Magnetogram 01.00(UT) GG1.88: 12.00(UT)GG18.88, Acta XII GNGTS, Roma.
 Di Gennaro, E., Baralli, F., Bovio, E., Faggioni, O. and Soldani, M. (2008) Clearance Operation of Teulada Site (Italy): A Novel Approach for Short Term MCM Missions in Sea Flor Hard Condition, Undersea Defence Technology UDT XXI, Glasgow.
 Faggioni, O., Soldani, M., Gabellone, A., Maggiani, P.V. and Leoncini, D. (2008) Development of Anti Intruders Underwater System: Time Demain Evaluation of the Self-Informed Magnetic Networks Performance. In: Corchado, E., Zunino, R., Gastaldo, P. and Herrero, á., Eds., Proceedings of the International Workshop on Computational Intelligence in Security for Information Systems CISIS’08. Advances in Soft Computing, Springer, Berlin, Heidelberg, 100-107.