Laurent Bouri

Dominique Lavenier

Evaluation of long read error correction software

Keywords: 3rd generation sequencing, long reads, correction séquençage de 3ème génération, lectures longues, correction

This report compares several read error correction software that use 3rd generation sequencing technology (long reads). The experimentations have been performed on several reference genomes and the results evaluated with QUAST. The long read error correctors that have been evaluated are: LSC-2, Proovread, Ectools, Lordec, Nanocorr, Nas, Jabba, Pacbiotoca, Lorma et MHAP. The first 8 software can merge long and short reads, while the last 3 software use only long reads.

1 Introduction

Contexte

Third-technology sequencing brought by Pacbioscience or Oxford Nanopore is able to produce average long read length of more than 10,000 bp and thus can advantageously be used to improve the genome assembly. In fact, long reads span more repetitive elements and thus can produce more contiguous reconstruction of the genome. However, long reads have raw error rate ranging from 10% to 15%, requiring a preliminary stage of correction before the assembly process.

The available correction software are mainly based on two strategies :

• hybrid correction (both long reads and short reads are required) ;

• denovo correction (only long reads are required) ; Hybrid correction uses short reads, such as Illumina, which have a much lower error rate, to correct long reads. The second strategy consists in aligning long reads against themselves. The evaluated hybrid correctors are listed in the table below. Some of them have the possibility of both recovering only the corrected regions ('trim' function) or the corrected and non-corrected regions ('untrim').

Evaluated long read error correctors

Correcteurs

trimmed reads untrimmed reads LSC 2 [al12a] yes yes Pacbiotoca [al12b] no yes Ectools [al14] no yes Proovread [T14] yes yes Lordec [START_REF] Rivals | LoRDEC: accurate and efficient long read error correction[END_REF] yes yes Nanocorr [al15c] no yes Nas [al15b] no yes Jabba [al16a] yes no 1.2.2 De-novo correctors (Long Read Only)

The following de-novo correctors have been evaluated:

• Pacbiotoca [al12b]

• Lorma [al16b] • MHAP (CANU) [al15a] 1.3 Method

Evaluation of long read error correction

In order to show the efficiency of long read error correction, several metrics reflecting long read quality can be calculated before and after assembly.

First, BLASR generates alignment of corrected long reads against the reference genome. From these alignments the following metrics are extracted:

• The average proportion of long read regions, aligned to the reference genome.

• The average number of match per long read and the percentage of identity of the aligned regions.

• The sum of corrected read length divided by the uncorrected read length. This gives an idea of the amount of sequences that has actually been corrected.

Then, the corrected long reads are assembled by Smartdenovo, a de-novo assembler giving good results, even for large genomes (> 100Mb) or with low coverage. Lastly, the QUAST (QUality ASsessment Tool) software [al13] evaluates the assembly by extracting several metrics:

• The number of contigs,

• Genome assembly length,

• N50,

• The fraction of the reference genome found among the contigs. This fraction is deduced by aligning the assembled genome against the reference genome using the MuMmer software [al99].

Datasets

The table below shows the 4 reference genomes and the datasets used for the various tests.

Input data

LSC-2 takes FASTA or FASTQ files as input.

Pipeline

runLSC script divides long read error correction into 5 steps:

• The sequences of long and the short reads are transformed by homopolymer compression so that each sequence of the same nucleotide is replaced by a single nucleotide of the same type.

• The short reads quality is checked. Indeed, some of these reads contain too much N letters or are too short.

• The short reads are aligned against the long ones with Bowtie2.

• The long reads are then modified according to the information provided by the short read consensus obtained with previous alignment.

• Once the correction points have been replaced by the corresponding short reads consensus, the rest of the compressed points are decompressed.

runLSC.py --long_reads LR.fa --short_reads SR.fa --specific_tempdir temp --output output_dir

• longs reads : long read file.

• short reads : short read file.

• specific tempdir : folder containing temporary files (optional).

• output : final assembly folder.

Inria

Encountered errors

Died at /LSC-2.0/bin/../utilities/explode fasta.pl ValueError: invalid literal for int() with base 10: " solution : Convert short reads from FASTQ format to FASTA format (AWK command line or Biopython Seq.IO module).

[bam header read] EOF marker is absent. The input is probably truncated.

solution : Install Bowtie2.

Output data

The corrected sequences are written into the "corrected read.fasta" file, while full LR.fasta file contains concatenate uncorrected terminus sequences and corrected sequences. Both files are located in the final assembly folder.

PacBioToCA

Introduction

Pacbiotoca is an hybrid and de-novo corrector, taking as input long reads in FASTQ file format and short reads in a frg file.

Website: http://wgs-assembler.sourceforge.net/wiki/index.php/PacBioToCA

Installation

PacBioToCA can be downloaded as source code or pre-compiled binaries.

Code source compilation:

Input data

PacBioToCA requires as input an estimation of the genome length, a configuration file, the long reads in a FASTQ file format and a fragments file "frg" containing short reads in the case of an hybrid correction. Programs such as fastatoCA and fastqtoCA are available in the PBcR package to perform conversions from common formats such as FASTA or FASTQ. Converting the FASTA file to .frg requires not only a file containing the sequences but also a quality value FASTA file.

fastaToCA -l libraryname -s seq.fasta -q qlt.fasta > seq.frg

• l : library name

• s : short reads in FASTA file format

• q : quality values FASTA file fastqToCA -libraryname LIB -technology pacbio-corrected -reads seq.fastq > seq.frg

• libraryname : library name

• technology : type of data (pacbio, illumina, 454,...)

• reads : short reads in FASTQ file format

Inria

Pipeline

PacBioToCA is able to perform hybrid or de-novo correction. In order to launch a de-novo error correction, simply do not provide the .frg file (short reads).

./pacBioToCA -threads 8 -libraryname name -s spec_file.spec -genomeSize <genome length> -fastq Long_reads.fastq short_reads.frg

• libraryname : prefix used for naming output files

• s : file containing specific options for correction

• genomeSize : genome length

• fastq : long reads to correct in FASTQ file format.

Output data

The corrected reads are stored in the two output files libraryname.fasta and libraryname.fastq.

Ectools

Input data

Ectools needs a contig file generated from short reads as well as a file including long reads.

Pipeline

Set the variable ECTOOLS_HOME with the path for input data:

$ ECTOOLS_HOME=/path/to/this/directory

Create a folder for long read error correction:

$ mkdir organism_correct

Create symbolic link to contig file:

$ cd organism_correct $ ln -s /path/to/organims_contigs.fasta

Ectools authors recommend keeping reads larger than 1Kb. Then, partition long reads:

$ python ${ECTOOLS_HOME}/partition.py 20 500 pbreads.length_filtered.fa

Copy the bash script for correction to the main folder:

$ cp ${ECTOOLS_HOME}/correct.sh .

Inria

Modify the global variables at the beginning of the correct.sh script and then run the next command line to correct long reads: $ for i in {0001..000N}; do cd $i; qsub -cwd -j y -t 1:M ../correct.sh; cd ..; done

• N : number of folders created by the "partition.py" script

• M : number of files in each folder.

Several steps are then performed to correct the long reads:

• Alignment of Pacbio reads against contigs with Nucmer.

• Selection of a set of contigs that covers best each long read to be corrected

• The "show-snps" program identifies the nucleotides that differ between contigs and long reads.

• Correction of long reads (Python script)

Output data

The corrected sequence files are written in the specified output folder: "p000N.cor.fa". Finally, merge files with a bash command:

$ cat p00*.cor.fa > organism.cor.fa

proovread

Introduction

Proovread is a de novo corrector, using a de-Bruijn graph constructed from long reads.

Website : https://github.com/BioInf-Wuerzburg/proovread

Installation

Proovread is available on Linux and requires NCBI Blast-2.2.24+, samtools-1.1+, Perl 5.10.1+ and the perl moduls Log::Log4perl and File::Which.

Code source compilation:

$ git clone --recursive https://github.com/BioInf-Wuerzburg/proovread $ cd proovread/util/bwa $ make

Input data

In order to not overload the processor and the memory, it is wise to divide the long reads data (FASTA, FASTQ) into several files:

$ SeqChunker -s 20M -o 0%03d pacbio_file
• s : file length

• o : output file name

Pipeline

Run the long read error correction process with the binary "proovread", located in the folder "bin", for every folders created in the previous step:

$ for i in {0001..000n}; do proovread -l $i -s /Path/to/short_reads {pre pb_$i; done

• n : Number of files generated by SeqChunker

• l : raw noisy long reads

• s : file containing accurate short reads

• pre : prefix used to name output file.

It is also possible to add unitigs with the argument "-unitigs".

Proovread corrects long reads in 2 steps:

Inria 1.
The mapping of short reads on long reads is done by SHRIMP2, by adapting the score mode to consider that insertions are more frequent than deletions and that substitutions are rare events. Bowtie2 and bwa mem are also supported.

2. A consensus sequence is computed from these alignments.

Output data

The corrected sequences are written in the specified output folder (trimmed and untrimmed reads)

Nanocorr

Introduction

Nanocorr is a hybrid corrector designed for nanopore long reads using blastn to align short reads on long reads.

Website: https://github.com/jgurtowski/nanocorr

Installation

Nanocorr has been designed to work in a SGE or similar environment (GNU parallel). Note that the python module "h5py" requires the installation of the hdf5 library Import nanocorr:

$

Input data

Nanocorr is supposed to correct only long reads from the nanopore technology (FASTA or FASTQ file format).

Pipeline

First, partition long reads into multiple files:

Inria $ python partition.py 100 500 nanopore_reads.fa

Then, run the "nanocorr.py" script in order to start long read error correction:

$ qsub -cwd -v PATH,LD_LIBRARY_PATH -t 1:500 -j 500 \ -o nanocorr_out /path/to/nanocorr.py query.fa reference.fa

• t : Number of files in the folder containing the partitioned long reads

• j : Declares whether the standard error stream of the task will be merged with the standard output stream of the same task

Nanocorr launches blast to align short reads on Nanopore long reads and then uses a dynamic programming algorithm based on the longest increasing subsequence problem to select sets of short reads corresponding to each long read. The consensus reads are then calculated using pbdagcon

Output data

The corrected sequences is divided into several files "p000N.blast6.r.fasta" and need to be concatenated in a single file:

LoRDEC

Introduction

LoRDEC is an hybrid corrector, using a de-Bruijn graph constructed from short reads to correct long reads.

Website : http://www.atgc-montpellier.fr/lordec/

Installation

LoRDEC is available on linux and requires Cmake 2.6+ and GCC 4.7+.

Import LoRDEC and the GATB library (http://gatb-core.gforge.inria.fr/) :

$ wget http://www.atgc-montpellier.fr/download/sources/lordec/LoRDEC-0.6.tar.gz $ tar zxvf LoRDEC-0.6.tar.gz $ cd LoRDEC-0.6 $ wget https://github.com/GATB/gatb-core/releases/download/v1.1.0/ \ gatb-core-1.1.0-bin-Linux.tar.gz $ tar zxvf gatb-core-1.1.0-bin-Linux.tar.gz Modify the variable GATB VER from the Makefile (1.1.0) Install LoRDEC $ make $ cd ..

Input data

LoRDEC requires short reads in FASTA or FASTQ file format and long reads in FASTA or FASTQ file format.

Pipeline

Run the long read error correction with the binary "lordec-correct":

$ lordec-correct -2 illumina.fasta -k 19 -s 3 -i pacbio.fasta \ -o pacbio-corrected.fasta
• 2 : File of short reads.

• k : Size of the kmer used in the de-Bruijn graph

• s : Abundance threshold of a kmer to be considered correct

• i : Input file • o : Output file Inria
A series of steps is then performed in order to correct the long reads:

1. Construction of a de-Bruijn graph from the short reads 2. Suppression of k-mer with occurrence less than the s value 3. Choose an optimal path of the graph by calculating the edit distance between the path and a region of long read.

Output data

The corrected sequences will be in the output file indicated after the "-o" parameter. The output file in FASTA format contains long reads. Corrected sequences are defined by uppercase letters while uncorrected sequences appears as lowercase letters. Lordec offers the possibility to remove the uncorrected sequences at the beginning and at the end of the long reading or to keep only the corrected sequences.

$ lordec-trim -i fichier_reads.fasta -o fichier_trim.fasta

• i : corrected reads file

• o : output file $ loredec-trim-split -i fichier_reads.fasta -o fichier_trim_split.fasta
• i : corrected reads file

• o : output file

Nas

Introduction

Nas adopts an hybrid approach to correct long reads from nanopore technology. The mapping of illumina sequences on nanopore reads will give long reads without errors called NaS (Nanopore Synthetic-long).

Website : https://github.com/institut-de-genomique/NaS

Installation

Nas is designed to run on linux or on grid environments (SGE,...). It can also be installed in a docker container (http://registry.hub.docker.com/u/rdbioseq/nas/). Nas requires the following dependencies:

• Shell tool GNU Parallel ( http://www.gnu.org/software/parallel/ ) at least (22052015)

• Perl 5.8.0 or higher

• Perl graph library (http://search.cpan.org/ jhi/Graph/)

• Perl GetOpt module (http://search.cpan.org/dist/Getopt-Long/)

• Perl Set::IntSpan module (http://search.cpan.org/ swmcd/Set-IntSpan-1.19/IntSpan.pm)

• Newbler assembler v2.9 (available from 454 website)

• Blat binary (at least v35) accessible through your PATH environnment variable

• Last binary (at least 588) accessible through your PATH environnment variable -Nas Installation:

$ git clone https://github.com/institut-de-genomique/NaS.git Specify the paths to Last and Blast in the "NaS" binary located in the "NaSv2" folder

Input data

Nas is able to correct reads smaller than 60kb and whose genome is smaller than 30Mb

Pipeline

Run long read error correction:

$ ./NaS --fq1 short_reads1.fastq --fq2 short_reads2.fastq --nano Longs_reads.fa \ --out NaS_output --mode sensitive --nb_proc 5

Inria

• nano : input file containing long reads

• mode : mapper used (fast=blat; sensitive=last)

• nb proc : number of processor Nas corrects long reads in two steps:

1. First, alignments are made between short reads and long reads (BLAT (fast mode) or LAST (sensitive mode)) in order to identify "seed-reads". The number of these "seed-reads" is then extended by looking for similar short reads in the initial set using comparisons.

2. Then, for each long reads, a microassembly is made from the short reads recruited using an Overlap-Layout-Consensus (Newbler) strategy.

Output data

The corrected sequences are written into the specified output folder, contained in the "NaS hqctg reads final.fa" file.

Jabba

Introduction

Jabba is an hybrid corrector, taking as input a de-Bruijn graph resulting from an assembly of short reads and uses it to correct long reads.

Website : https://github.com/gmiclotte/Jabba

Installation

Jabba is available for Linux and requires Cmake 2.6+ et GCC 4.7+. It is possible to install Karect for the step of short read correction and Brownie for the construction of the de-Bruijn graph. Note that Jabba needs sparsehash to be installed.

Install karect:

$ git clone https://github.com/aminallam/karect.git $ cd karect $ make $ cd ..

Install sparsehash

$ git clone https://github.com/sparsehash/sparsehash.git $ ./configure {prefix=<chemin vers le dossiersparsehas> $ make

Install brownie:

$

Installer Jabba:

$ git clone https://github.com/gmiclotte/Jabba.git $ cd Jabba $ ./compile.sh $ cd ..

Inria

Input data

Jabba takes long reads in a FASTA or FASTQ file format and the de-Bruijn graph constructed by Brownie from the short reads, themselves corrected by a tool like Karect.

Pipeline

First, the authors recommend to start the process with karect:

$ mkdir karect_output $ ./karect/karect -correct -matchtype=hamming -celltype=haploid \ -inputfile=short_reads.fastq -resultdir=karect_output -tempdir=karect_output

Then run long read error correction as follow:

$ /Jabba/bin/jabba -o jabba_output -k 75 -brownie_output/DBGraph.fasta \ -fastq long_reads.fastq

• l : minimum seed length

• k : kmer length used in Debruijn graph.

However, for each new tested dataset and the construction of a new de-Bruijn graph, the corrected reads must be redirected to a new output file. Jabba corrects the long reads by aligning them on a de-Bruijn graph using a "seed and extend" method where the "maximum exact matches" represent the seeds.

Output dataset

The corrected sequences will be in the specified output folder, contained in a "Jabba-input filename.fasta" file.

3 Denovo correctors (Long Read only)

LoRMA

Introduction

LoRMA is a de novo corrector using a de-Bruijn graph constructed from long reads to correct themselves.

Website : https://www.cs.helsinki.fi/u/lmsalmel/LoRMA/

Installation

LoRMA is available on Linux and requires Cmake 2.6+, GCC 4.7+ and Lordec.

Import LoRDEC and the GATB library:

$ wget http://www.atgc-montpellier.fr/download/sources/lordec/LoRDEC-0.6.tar.gz $ tar zxvf LoRDEC-0.6.tar.gz $ cd LoRDEC-0.6 $ wget https://github.com/GATB/gatb-core/releases/download/v1.1.0/ \ gatb-core-1.1.0-bin-Linux.tar.gz $ tar zxvf gatb-core-1.1.0-bin-Linux.tar.gz Modify the GATB VER variable from Makefile (1.1.0 here) Install LoRDEC $ make $ cd ..

Install loRMA

$ wget https://www.cs.helsinki.fi/u/lmsalmel/LoRMA/LoRMA-0.3.tar.gz $ tar zxvf LoRMA-0.3.tar.gz $ cd LoRMA-0.3 $ mkdir build $ cd build $ cmake .. $ make Change the LORDECDIR variable in the lorma.sh file to specify the path to the LoRDEC installation folder.

Input data

LoRMA takes long reads in a FASTA or FASTQ file format as input.

Inria

Pipeline

LoRMA software correct long reads in two steps:

1. Iterative construction of a Dde-Bruijn graph from long reads, with an increasing value of the kmer size. Each node in the graph represents a genomic kmer, selected according to its abundance in the dataset of long reads. For each long read, identify possible paths in the graph between the solid kmer and then select the best path. Correct the regions of the long read.

2. Finally, use the graph to identify similar reads. Align these reads in order to build a consensus.

Output data

The corrected sequences are located in a file named "final.fasta"

MHAP (CANU)

Introduction CANU includes a de novo correction module, using MHAP to align long reads against themselves and pbdagcon to correct the long reads by a consensus step.

Website : https://github.com/marbl/canu

Installation

Canu is designed to run on grid environments (LSF/PBS/Torque/Slrum/SGE are supported) and requires java 1.8 and gnuplot.

Install Canu :

$ git clone https://github.com/marbl/canu.git $ cd canu/src $ make -j <nombre de threads>

Input data

Canu Requires only long reads in FASTA or FASTQ format.

Pipeline

Run the long read error correction with the "canu" binary, located in the folder "canu/Linuxamd64/bin":

$ ./canu -correct -p ecoli -d ecoli genomeSize=4.8m -pacbio-raw longs-reads.fastq

• genomeSize : Genome size

• pacbio-raw : Type of long reads (pacbio-raw or naopore-raw)

It is also possible to isolate the corrected sequences with the "canu" binay:

$ ./canu -trim -p ecoli -d ecoli genomeSize=4.8m -pacbio-corrected corrected-reads.fastq

Canu corrects long reads in two steps:

1. At first, MHAP finds inaccurate overlaps between long reads.

2. Pbdagcon then builds a consensus sequence from these overlaps.

Encountered errors

canu failed with 'can't configure for SGE' solution : Add the following arguments in the command line: gridEngineMemoryOption="-l h vmem=MEMORY" gridEngineThreadsOption="-pe make THREADS"

Inria

Output data

The corrected sequences will be stored in multiple FASTA files concatenated in the output folder: "correction/2-correction/correction outputs"

4 Evaluation of correction and assembly

Evaluation methods

Long reads correction will be evaluated before and after assembly. First, a series of information is extracted either directly from the corrected long reads or from alignments performed with BLASR between long reads and the reference genome:

• Size : The sum of the length of the corrected long reads, divided by the length of uncorrected long reads.

• Execution time : Corrector execution time.

• % mapped region : Proportion of corrected reads that BLASR succeeded in mapping to the reference genome.

• Mean number of match : The average number of match per long read.

• % identity : Percentage of identity given by BLASR for the region of long read that could be aligned.

Then, the reference genomes allow to evaluate the quality of the different assemblies, based on the results of the metrics produced by the QUAST software. The different metrics are listed below:

• # contigs (> 1000pb) : Total number of contigs exceeding 1000 bp after assembly.

• Total length : Assembly length.

• N50 : Length for which the collection of all contigs of this length or greater, covers at least half of the assembly.

• Genome fraction (%) : the percentage of aligned bases in the reference genome. The lack of correction by the Nas tool is due to an incompatibility with the illumina data. The lack of correction by the Nas tool is due to an incompatibility with the illumina data.

de-novo correctors

de-novo correctors

Inria

The lack of results after the assembly of long reads corrected by loRMA is due to a low quantity of available long reads after the correction stage. The lack of correction of the Nas tool is due to an incompatibility of the illumina data. The lack of results after the assembly of long reads corrected by PacBioToCA ou Jabba is due to a low quantity of available long reads after the correction stage.

de-novo correctors

Inria

1. 2 . 1

 21 Hybrid correctors (Short and Long Reads)

$

  cat p*.blast6.r.fasta > output_file.fasta Encountered errors c.fatal error: hdf5.h: No such file or directory virtualenv solution : installation of hdf5 library

•

  Oxford NanoPort reads (MinIon) corrected by Lordec (10x coverage) : 3427 reads • Illumina reads (MiSeq) : 211219 reads of length 2x150pb

4. 9 •

 9 Test 6 : Escherichia coli k-12, reads Pacbio 100x (P6-C4) Dataset: Pacbio reads corrected by lordec (100x coverage) : 87497 reads • Illumina reads (MiSeq) : 16759877 reads of length 2x150pb • Contigs generated by sparse assembler from Illumina reads ( 1876792 contigs)

  Since most of the long read error correctors requires a large amount of hardware resources (CPU, RAM), their evaluation was done by submitting tasks to the GenOuest 1 platform cluster.

	1.3.3 Hardware resources
	Correctors	Number of threads
	LSC 2	8
	PacbiotoCa	8
	Ectools	1
	Proovread	8
	Lordec	8
	Nanocorr	1
	Nas	8
	Jabba	8
	MHAP	8
	Lorma	8
	Cluster node configuration:
	Genome	length (M bp) Test	Minion	PacBio	Illumina
	Acinetobacter 3.6 M	1	10x, 3.4K reads	211K reads
				2	20x, 10K reads	211K reads
	E.Coli	4.6 M	3		10x, P4C2, 36K reads	16M reads
				4		100x, P4C2, 91K reads	16M reads
				5		10x, P6C4, 8.7K reads	16M reads
				6		100x, P6C4, 87K reads	16M reads
				7	20x, 22K reads	16M reads
	S. Cerevisa	11.6 M	8		10x, P4C2, 26K reads	3.8M reads
				9		100x,P4C2, 261K reads	3.8M reads
				10	20x, 47K reads	3.8M reads
	C. Elegans	100 M	11		10x, P6C4, 92K reads	55M reads
				12		100x, P6C4, 740K reads	55M reads
						Inria

• Number of CPU: 40 • CPU frequency: 2.6 GHz • Available RAM: 256 GBytes 1 GenOuest: https://www.genouest.org/ RR n°9028

  IntroductionLSC-2 is an hybrid corrector of long reads. Long reads and short reads are first compressed into homopolymers, then short reads are mapped to long reads with Bowtie2. Finally, the short read consensus replaces the long read sequences.

	2 Hybrid correctors (short and long reads)
	2.1 LSC-2

Website: http://www.healthcare.uiowa.edu/labs/au/LSC/default.asp Installation LSC-2 can be downloaded as pre-compiled binaries. Bowtie2 is required and must be installed.

Extraction of pre-compiled binaries files:

$ tar zxvf LSC-2.0.tar.gz

•

  Contigs generated by sparse assembler from Illumina reads (233 contigs)

	4.5 Test 2 : Acinetobacter sp, ADP1, run6; Minion 20x 4.6 Test 3 : Escherichia coli k-12, reads Pacbio 10x (P4-C2)
	Dataset: Dataset:
	• Pacbio reads corrected by lordec (10x coverage) : 36355 reads • Oxford NanoPort reads (MinIon), corrected by lordec (20x coverage) : 10116 reads • Illumina reads (MiSeq) : 16759877 reads of length 2x150pb
	• Illumina reads (MiSeq) : 211219 reads of length length 2x150pb • Contigs generated by sparse assembler from illumina reads (1876792 contigs)
	• Contigs generated by sparse assembler from Illumina reads (233 contigs)
	Hybrid correctors Hybrid correctors hybrid correctors
	LSC-2 untrim 88.04 5 342s 94.76 10 135 trim 79.14 5 342s 99.11 10 731 87.49 86.68 29 34 153 908 120 006 3 267 123 2 668 193 15.19 13.62 LSC-2 PacBioToCA PacBioToCA 16.42 1 084s 99.25 1 904 93.07 7 19 707 153 749 0 Ectools LSC-2 PacBioToCA untrim trim untrim 90.39 62.25 81.30 86.11 95.40 80.38 8 935s % mapped region Metrics Size Execution time Mean number of match % identity contigs (>=1000pb) N50 Total length Genome fraction Metrics trim Metrics Size 81.75 Size Execution time 8 935s 4 088 120s Execution time 31 768s 31 768s 17 298s % mapped region 99.15 95.45 99.63 99.7 % mapped region 94.93 90.23 98.01 Mean number of match 10 005 10 561 1 227 12 436 Mean number of match 4 266 4 433 1 596 % identity 88.73 87.88 99.93 99.6 % identity 95.47 94.35 99.85 contigs (>=1000pb) 1 2 47 1 contigs (>=1000pb) 71 66 40 N50 3 605 124 2 415 855 85 022 3 649 197 3 560 442 879 221 3 446 093 Ectools Proovread trim untrim 27.04 40.39 101.02 90s 5 214s 5214s 99.7 92.69 34 12 417 6 442 11 239 99.6 99.93 98.37 17 20 18 338 569 50 752 315 675 96.75 24.42 92.54 Ectools Proovread Proovread trim trim untrim untrim 32.07 80.39 79.18 106.9 83.84 7 087s 160s 2 791s 2 791s 7 087s 99.8 99.6 98.51 91.04 97.05 4 559 6 188 3 634 4 631 11 831 99.95 99.3 99.82 97.64 98.99 29 61 65 77 1 175 774 N50 41 815 42 879 25 658 35 761 32 061 44 775 3 631 367 Total length 3 605 124 3 598 621 2 630 061 Total length 2 723 468 2 526 333 1 000 670 1 935 279 2 069 936 3 005 248 3 649 197 3 063 864 3 631 367 Genome fraction 25.82 28.88 73.07 98.79 84.98 97.21 Genome fraction 17.56 50.05 21.55 41.52 44.40 63.44
	LoRDEC trim untrim 79.91 101.82 521s 521s 99.81 91 6 607 11 038 99.9 97.13 29 26 178 054 164 126 3 369 809 3 341 599 3 329 760 3 203 245 2 809 569 Nanocorr Nas Jabba 72.53 96.1 82.9 1 749s 22 500s 24.72s 99.33 99.45 88.8 9 668 13 320 8 658 97.06 99.93 99.9 27 20 39 127 495 211 411 90 115 92.87 89.02 87.31 88.76 78.08 LoRDEC Nanocorr Nas LoRDEC Nanocorr Nas Jabba Jabba trim Mean number of match Metrics Size Execution time % mapped region % identity contigs (>=1000pb) N50 Total length Genome fraction Metrics Metrics trim untrim untrim Size 82.28 102.02 76.14 98.33 Size 73.24 94.67 80.95 -80 84.4 Execution time 661s 661s 2 520s 91 920s Execution time 2 392s 2 392s 15 480s -267s 29s % mapped region 99.82 92.51 99.32 99.32 % mapped region 99.10 93.16 96.53 -99.9 99.8 Mean number of match 6 575 10 919 9 813 12 808 Mean number of match 1 034 2 517 4 150 -4 803 8 658 % identity 99.92 97.20 97.39 99.94 % identity 99.53 94.82 97.95 -99.9 99.9 contigs (>=1000pb) 4 1 1 32 contigs (>=1000pb) 1 90 67 -77 28 N50 1 113 078 3 622 997 3 647 606 116 593 N50 12 447 33 149 44 366 -46 908 110 934 Total length Total length 12 447 2 787 937 2 573 103 -3 155 083 3 601 490 3 622 997 3 647 606 3 230 285 2 648 543 Genome fraction 99.99 97.22 94 89.69 73.59 Genome fraction 0.26 53.04 53.78 -68
	de-novo correctors de-novo correctors
	Metrics	Metrics	PacBioToCA MHAP (CANU) LoRMA PacBioToCA MHAP (CANU) LoRMA
	Size	Size			68.47	30.44	54.01	33.36	21.1	0.37
	Execution time Execution time	6 745s	1 037	540s	120	670s	243.84
	% mapped region % mapped region		99.15	99.25	99.67	99.6	99.6	99.6
	Mean number of match Mean number of match 7 297	1 869	8 631	6 381	189	174
	% identity % identity			93.10	93.10	99.49	98.83	98.13	98.58
	contigs (>=1000pb) contigs (>=1000pb)	14	6	16	23	no data	no data
	N50	N50	2 020 356	23 350	1 892 568	29 225 no data	no data
	Total length Total length	3 249 142	132 365	3 230 123	628 069 no data	no data
	Genome fraction Genome fraction		1.3	0	7.93	0	no data	no data
									Inria

  The lack of correction by the Nas tool is due to an incompatibility with the illumina data.

	4.8 Test 5 : Escherichia coli k-12, reads Pacbio 10x (P6-C4)
	Dataset:
	• Pacbio reads corrected by lordec (couverture 10x) : 8746 reads
	• Illumina reads (MiSeq) : 16759877 reads of length 2x150pb
	• Contigs generated by sparse assembler from Illumina reads(1876792 contigs)
	hybrid correctors
					LSC-2		PacBioToCA	Ectools	Proovread
		Metrics			trim	untrim		trim	untrim
		Size			91.70	94.60		87.2	77.88	75.06	79.32
	Execution time	110 418s	110 418s	45 542s	60s	12 840s	12 840s
	% mapped region		97.23	95.03		98.34	99.53	99.34	93.56
	Mean number of match		7 471	7569		2 756	10 144	7 054	9 121
		% identity			97.30	96.93		99.67	99.2	99.85	98.96
	contigs (>=1000pb)		6	7		-	20	7	7
		N50		1491698	1 308 005	No data	525 979	1 486 401 1 259 739
	Total length		4 676 830 4 677 980	No data	4 478 726 4 639 996 4 663 929
	Genome fraction		97.15	96.2		No data	95.98	99.38	99.17
					LoRDEC		Nanocorr Nas	Jabba
		Metrics			trim	untrim
		Size			84.88	96.62		67.31	-	86.4
	Execution time		603s	603s		20 460s	-	274s
	% mapped region		98.84	93.64		96.87	-	99.8
	Mean number of match		1 976	7 534		5 942	-	7 466
		% identity			99.58	97.27		98.37	-	99.99
	contigs (>=1000pb)		79	9		34	-	48
		N50			31 940	1 207 236	215 605	-	99 142
	Total length		2 323 363 4 650 012 4 443 229	-	4 103 224
	Genome fraction		50.03	94.86		92.40	-	88.4
	de-novo correctors
	Metrics	Metrics	PacBioToCA MHAP (CANU) PacBioToCA MHAP (CANU) LoRMA LoRMA
	Size Execution time Size Execution time % mapped region % mapped region Mean number of match Mean number of match 38.63 3 219s 99.19 5 994 % identity 98.47 % identity	76.04 2 738s 99.53 5 336 99.18	29.55 3 180s 96.01 7 172 99.84	67.9 2 340s 99.3 7 396 99.3	79.6 7 524s 99.7 1 609 99.5	5.3 323s 99.7 283 99.68
	contigs (>=1000pb) contigs (>=1000pb) N50 4 624 650 1 N50 Total length 4 624 650 Total length 3 318 200 57 66 797 4 647 081 1 22 729 3 710 724 50 100 017 86 784 70 4 109 490 Genome fraction 99.96 9.28 88.3 Genome fraction 71.62 80.13	----

  The lack of correction of the Nas tool is due to an incompatibility of the illumina data.The lack of correction of the Nas tool is due to an incompatibility of the illumina data. The lack of results after the assembly of long reads corrected by PacBioToCA ou Jabba is due to a low quantity of available long reads after the correction stage.The lack of correction of the Nas tool is due to an incompatibility of the illumina data.Inria4.13 Test 10 : Saccharomyces cerevisiae W303, reads Minion 20xThe lack of correction of the Nas tool is due to an incompatibility of the illumina data. The lack of results after the assembly of long reads corrected by PacBioToCA ou Jabba is due to a low quantity of available long reads after the correction stage.

	4.10 Test 7 : Escherichia coli k-12, reads Minion 20x 4.11 Test 8 : Saccharomyces cerevisiae W303, reads Pacbio 10x (P4-4.12 Test 9 : Saccharomyces cerevisiae W303, reads Pacbio 100x (P4-
		C2) C2)
	Dataset: The lack of results after the assembly of long reads corrected by loRMA is due to a low quantity of Dataset:
	• Minion reads corrected by lordec (20x coverage) : 22270 reads Dataset: available long reads after the correction stage. Dataset:
	• Illumina reads (MiSeq) : 16759877 reads of length 2x150pb • Pacbio reads corrected by lordec (100x coverage) : 261964 reads • Minion reads corrected by lordec (20x coverage) : 47027 reads
	• Contigs generated by sparse assembler from Illumina reads( 1876792 contigs) • Pacbio reads corrected by Lordec (10x coverage) : 26196 reads • Illumina reads (HiSeq) : 3815678 reads of length 2x100pb • Illumina reads (HiSeq) : 3815678 reads of length 2x100pb
	• Illumina reads (HiSeq) : 3815678 reads of length 2x100pb • Contigs generated by sparse assembler from Illumina reads (10055 contigs) • Contigs generated by sparse assembler from Illumina reads (10055 contigs)
	hybrid correctors • Contigs generated by sparse assembler from Illumina reads (10055 contigs)
	hybrid correctors		LSC-2		PacBioToCA	Ectools	Proovread
	Metrics			trim	untrim				trim	untrim
	Size hybrid correctors hybrid correctors	88.49	99.87 LSC-2	90.34	80.39 PacBioToCA	94.3 Ectools	100.9 Proovread
	Execution time Metrics	35 918s	35 918s trim	40 138s untrim	160s	3 365s	22 140s trim	untrim
	% mapped region Mean number of match Metrics Size Execution time Metrics	99.16 5 428 LSC-2 trim 95.76 5933 LSC-2 untrim 98.39 1 749 PacBioToCA 99.9 7 838 83.24 93.84 15.53 PacBioToCA Ectools 943 079s 943 079s 33 179s trim untrim	99.78 4 277 Ectools 55.70 Proovread 95.47 5 881 Proovread trim untrim 49.56 96.94 1 325s 200 040s 200 040s trim untrim
	% identity % mapped region 88.35 Size Size 36.32	86.86 40.46 85.33 55.86	99.79 80.45 92.41 4.60		99.8 97.46 21.49 32.75	99.91 51.98 95.35 23.81	97.96 44.23 85.97 100.85	99.45 79.78
	contigs (>=1000pb) Execution time Mean number of match 5 Execution time 13 686s	6 77 167s 4 287 13 686s	36 10 866s 77 167s 4531		25 120s 11 172 802	2 6 719 147s 47 020s	1 28 260s 2 402 47 020s	28 260s 4 441
	N50 % mapped region 1 313 431 1 083 813 83.06 % identity 93.74 % mapped region 92.23 85.42	172 755 76.99 92.82 98.74	210 664 97.83 99.64 96.4	4 526 234 4 671 320 95.5 86.74 99.47 99.40 87.43 78.46	91.04 96.07
	Total length Genome fraction Mean number of match 4 645 507 4 699 525 17.56 17.29 4 154 % identity 92.45 contigs (>=1000pb) 23 Mean number of match 4 225 5921 N50 780551 % identity 79.67 77.40	4 510 878 97.20 4533 91.42 21 216 781 007 -	4 474 215 4 661 072 4 671 320 94.97 99.97 96.31 639 6 904 2 187 99.6 99.4 99.40 236 60 71 7 537 1 349 5 943 40 462 340 500 135 059 99.3 99.40 90.12	4 631 97.64 818 767 20
	LoRDEC 175 1 2182 386 12 181 497 Nanocorr Nas 169 -8 340 302 15 -140 Jabba trim contigs (>=1000pb) Total length contigs (>=1000pb) 76 Metrics N50 52 123 48 098 -51 046 173 11 628 184 7 019 792 12 178 144 79 173 6 131 25 102 66 148 Genome fraction 89.69 86.40 57.66 96.74 57.66 95.05 N50 25 411 23 975 -100 030 16 107 114 531 untrim Total length 7 677 266 7 234 422 -7 828 809 1 908 619 8 766 345 Total length 18 341 73 336 512 -10 456 921 91 934 11 046 053
	Size Genome fraction Genome fraction	87.54 0.15	101.08 19.81 LoRDEC 47.53 92.60 0.034 -	--Nanocorr 88.48 94.8 Nas Jabba 66.49 0.55	16 80.16	69.16
	Execution time Metrics		406s	406s trim	34 200s untrim		-	272s
	% mapped region Mean number of match Metrics Size Execution time Metrics	99.26 1 719 LoRDEC 94.67 5 763 LoRDEC trim untrim 98.46 5 505 76.77 97.66 Nanocorr 31 024s 31 024s trim untrim	--Nanocorr Nas Jabba 99.9 5 699 87.31 -Nas Jabba 11.6 143 100s -213s
	% identity % mapped region 99.63 Size Size 40.45	95.14 68.65 90.04 101.31	96.49 99.91 78.84 56.14		-42.41 90.35 -	99.9 -4.7	-	9.2 98.9
	contigs (>=1000pb) Execution time Mean number of match 114 Execution time 2 877s	2 1 546s 1 869 2 877s	1 1 546s 4 413 18 450s	-8 932s 3 227 -	51 -140s	-	122s 576
	N50 % mapped region 32 904 % identity % mapped region 94.64	4 650 531 4 673 315 90.91 76 99.41 96.39 78.74 95.96		-91.23 97.54 -	105 676 --99.6	99.1 99.9
	Total length Genome fraction Mean number of match 3 246 379 4 707 245 4 673 315 69.91 91.15 93.60 1 830 4 497 % identity 99.44 95.72 contigs (>=1000pb) 183 20 Mean number of match 1 269 6 062 6 542 N50 87 654 820 758 % identity 99.64 90.36 99.22		--2 842 97.25 24 -755 939 4 423 700 95.3 ---579 --99.9	574 99.9 --
	contigs (>=1000pb) Total length contigs (>=1000pb) 15	87 11 479 142 12 209 277 12 281 380 174 70 156 81 -	-	-	-	-	-
	N50 Genome fraction N50 16 772	20 734 92.70 84 777	60 876 95.26 200 230	24 708 97.32 -	-	-	-	-	-
	Total length Total length	1 820 797 11 633 571 1 744 892 244 645 10 456 229 11 362 908 -	-	-	-
	de-novo correctors Genome fraction Genome fraction	1.85		15.07 77.93	66.95 95.74		13.89 -	-	-	-
	Metrics de-novo correctors PacBioToCA MHAP (CANU) LoRMA
	Metrics Size Execution time Size Metrics Execution time % mapped region Size % mapped region Mean number of match 78.96 PacBioToCA MHAP (CANU) 76.99 13.4 30.20 21.16 6 396s 720s PacBioToCA MHAP (CANU) 238s 5 120s 7 380s 99.21 99.69 99.6 29.21 23.9 97.68 99.6 12 387 6 381 Mean number of match 4 923 5 577 164 % identity 92.76 99.24 98.53 de-novo correctors Execution time 10 916s 9 360s % mapped region 91.49 86.03 de-novo correctors % identity 98.60 98.8 contigs (>=1000pb) 9 1 -Mean number of match 6 923 8 327 contigs (>=1000pb) 1 1 N50 642 915 4 396 819 -Metrics PacBioToCA MHAP (CANU) LoRMA LoRMA LoRMA 83.9 70.3 18 923s 99.8 4 1036s 94.4 1 351 1 461 99.71 % identity 97.83 99.14 Metrics PacBioToCA MHAP (CANU) LoRMA 99.5 27 N50 4 646 758 4 646 370 Total length 4 321 437 4 396 819 -Size 16.87 13.93 6.17 contigs (>=1000pb) 23 24 210 Size 5.02 --306 808 Total length 4 646 758 4 646 370 Genome fraction 0.26 0.32 -Execution time 3 039s 6 660s 849s N50 771921 776764 64372 Execution time 4 908s --4 597 457 Genome fraction 99.97 99.99 The lack of results after the assembly of long reads corrected by loRMA is due to a low quantity of % mapped region 93.54 81.14 60.4 Total length 12 143 815 12 284 076 9 946 745 % mapped region 89.78 --98.9 available long reads after the correction stage. Mean number of match 1562 3 553 202 Genome fraction 97.99 98.53 84.9 Mean number of match 986 --
	% identity % identity		88.58	92.76	-	97.9		-	93.17
	contigs (>=1000pb) contigs (>=1000pb)	1	6	-	7		-	no data
	N50	N50		12991	17 647	-	24 044	-	no data	Inria Inria
	Total length Total length		12 991	91 430	-	158 844	-	no data
	Genome fraction Genome fraction	0.003	0.74	-	1.21		-	no data

RR n°9028

Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

The lack of results after the assembly of long reads corrected by PacBioToCA is due to a low quantity of available long reads after the correction stage.

Inria


Amar Zeher
email: amar.zeher@u-bourgogne.fr

Stéphane Binczak

Jerome Joli
email: jjoli@selecom.fr

Jérome Joli

Characterization of Radio Frequency Echo Using Frequency Sweeping and Power Analysis

Keywords: Radio frequency echo, gain margin, narrow band, frequency hopping

come

I. INTRODUCTION

In regions where reception quality is poor, it is usual to locate repeaters which are deployed to expand and enhance coverage region of Base Station Transceiver (BTS) costeffectively, without the need for modulation and demodulation setups. A repeater receives signals from the main transmitter which is generally BTS on the selected carriers ; then, it retransmits them to receivers which are generally Mobile Station (MS) such as cell phones, Personal Digital Assistant (PDA) and communication terminals. The receiving signal is filtered, amplified and re-transmitted; if the repeater retransmits signals in the same frequencies, it is called Iso-Frequency Repeater (IFR), otherwise, it is called Transposing Frequency Repeater (TFR). Deploying TFR to extend coverage area is robust, but it is expensive. Contrariwise, IFR are cost-effectively and simple to set up. Nevertheless, with IFR, if the channel insulation is insufficient, a part of the transmitted signal could be picked up by the receiving antenna of the repeater, leading to RFE creation.This phenomenon is well known as the issue was first posed in acoustics field [START_REF] Sondhi | Silencing echoes on the telephone network[END_REF]. When using repeaters, Echo Cancellation is one method of compensating for external Radio Frequency (RF) coupling between antennas [START_REF] Braithwaite | Echo cancellation for a wide bandwidth mixed-mode wcdma/gsm repeater with digital sub-band filtering[END_REF]. The echo signals must be reduced to an allowable level in order to avoid problems with distortion and oscillation [START_REF] Hong | Echo canceller for on-channel repeaters in t-dmb system[END_REF]. Many techniques exist for cancelling echoes, notably in the area of acoustics as Wiener Filter, Linear Prediction [START_REF] Hänsler | Acoustic Echo and Noise Control: A Practical Approach[END_REF] or Non-linear Acoustic Echo Cancellation Based on Volterra Filters [START_REF] Guerin | Nonlinear acoustic echo cancellation based on volterra filters[END_REF]. Some of these techniques are applied to radiocommunication. Depending on the nature of concerning signal, the efficiency of these methods could be called into question. Indeed, modulation constraints as guard interval, Error Vector Magnitude (EVM), Modulation Error Ratio (MER) have to be satisfied when filtering.

Narrow band signal (NBS) with frequency hopping, as in GSM or TETRA standards, is one major difficulty when using linear filtering because of the convergence rate of adaptive algorithms which are generally too slow when compared with the hopping frequency rate and randomness [START_REF] Lagrange | Réseaux GSM Des principes à la norme[END_REF]. Under these constraints, linear filtering increases transmission errors and creates an Adjacent-Channel Interference (ACI).

Radiofrequency signal can also be Wide Band Signal (WBS), which occupies a range of frequencies such as DVB-T/H (8 MHz in France) and CDMA. NBS has a narrower band width, generally less than 1 MHz as GSM (200 KHz) and TETRA standards (25 KHz). Part (a) of Fig. 1 shows an example of WBS and NBS displayed in the frequency domain without RFE.

Ripples indicates the presence of a RFE and, as illustrated in part (b) of the Fig. 1, they are only distinguishable in WBS and the noise floor parts. With noise floor or wide band, it is easily detectable and voidable by using linear adaptive algorithms.

Nevertheless, the narrower the band is, the more difficult these solutions can be used effectively. For example, adaptive filtering using Wiener solution based on the correlation calculation is a good solution when used with DVB-T or CDMA standard, but it cannot be used when GSM signal with frequency hopping is used, because of the ACI and the unwanted spurious which are created.

The transmitted signal can reach the receiving antenna by one or more paths (it is then called multipath), depending on the topology of the location. In this paper, RFE phenomenon is modeled using only the direct path. It is also supposed in the simulation part that the repeater includes three functions which are a controllable amplification, delay and Additive Gaussian White Noise (AWGN). The channel parameters are chosen to be a constant delay and variable attenuation (other parameters as fading are therefore not taken into account) and the propagation delay in wires is supposed negligible. The article is organized as follows: In section II, mathematical equations modeling the system are presented and the relationships between the power in the repeater and the system parameters as the repeater gain, repeater delay, channel attention are stated. These theoretical results are consistent with findings from the simulation presented in section III. Simulink tool is used to realize the closed-loop system and results are analyzed using Matlab. Finally, as a conclusion, the main results are summarized and discussed.

II. THEORY

Taking into account the set of assumptions previously mentioned, a system composed of a repeater, two antennas and a channel is simplified, as shown in Fig. 2. The closed Fig. 2. Radiofrequency coupling, considering only the direct echo path. loop of the system is expressed by the following equations :

s(t) = i(t) + e(t) , e(t) = ay(t -τ 1 ) , y(t) = gs(t -τ 0 ) , (1)
where i(t) is the incident signal, y(t), the transmitted signal and e(t), the feedback signal (echo). g and τ 0 are respectively the repeater gain and the repeater delay, while a and τ 1 are respectively the channel attenuation and the channel delay.

System (1) can be expressed such as

y(t) = g • i(t -τ 0 ) + ga • y(t -τ 0 -τ 1 ).
(2)

Eq. ( 2) is a recursive equation, which can be expressed by y(t), y(t -(τ 0 + τ 1 )), s(t -2(τ 0 + τ 1 )), . . . , s(t -l(τ 0τ 1 )), s(t -(l + 1)(τ 0 -τ 1 )), s(t -(l + 2)(τ 0 -τ 1 )), while l is the number of secondary echoes. Introducing

g • a = δ, τ 0 + τ 1 = τ,
eq. ( 2) can be expressed by the following sum

y(t) = gi(t -τ 0 )+ g l+1 n=1 δ n i(t -nτ -τ 0 + δ l+2 y(t -(l + 2)τ ) , (3)
which can be reduced to its simplest expression, so that

y(t) = g l+1 n=0 δ n i t -(nτ + τ 0 ) + δ l+2 y t -(l + 2)τ . (4)
Using Fourier Transformation, eq. ( 4) leads to

Y(ω) = δ l+2 e -jω(l+2)τ Y(ω) + gI(ω) l+1 n=0
δ n e -jω(nτ +τ0) , [START_REF] Guerin | Nonlinear acoustic echo cancellation based on volterra filters[END_REF] where Y(ω) and I(ω) are the Fourier Transformation of y(t) and i(t), respectively. The transfert function is defined as

H(ω) = Y(ω) I(ω) . ( 6
)
From eq. ( 5), H(ω) can be written such as

H(ω) = g l+1 n=0
δ n e -jω(nτ +τ0)

1 -δ l+2 e -jωτ (l+2) .

When the repeater steady state is reached, l tends towards infinity. Let H ∞ be this limit. Then, from equation (7), the system transfer function leads to

H ∞ (ω) = g e jω (τ +τ 0) -δ + e jω τ . (8)
Thus, the power p of the system can be written so that

p(ω) = |H ∞ (ω)| 2 = g 2 δ 2 -2 δ cos (ω τ ) + 1 . ( 9
)
Let us define in dB the following parameters : Using these notations with eq. ( 9), the power in dB is given by

P (ω) = G -10 log δ 2 -2δ cos (ω τ ) + 1 . (11)
In practice, most of the variables are expressed in dB, thus

P (ω) = G -10 log 10 ∆/10 -2 • 10 ∆/20 cos (ω τ ) + 1 . (12)
It is noted that all variables are real such that G ∈ R + , A ∈ R -, ∆ ∈ R. Note that in the case of strong echo, corresponding to the fact that gain margin tends towards zero or becomes positive, eq. ( 12) tends asymptotically to

P (ω) ≈ G -3 -10 log 1 -cos (ω τ ) . ( 13
)
Therefore, power value P can theoretically tend to infinity.

In the case of weak echo, corresponding to the fact that gain margin is negative and tends towards minus infinity, eq. ( 12) tends asymptotically to

P (ω) ≈ G. ( 14
)
It means that power is constant when sweeping frequency, as illustrated in Fig. 3. Now, let R be the magnitude of these ripples, defined as the difference between the maximum of power (P max ) and minimum of power (P min ). From eq (12), P max and P min can be written such as

P max (∆) = G -10 log 10 ∆/10 -2 • 10 ∆/20 + 1 , P min (∆) = G -10 log 10 ∆/10 + 2 • 10 ∆/20 + 1 , (15) yielding R(∆) = 20 log 10 ∆/20 + 1 10 ∆/20 -1 , (16)
which is indeed only depending on gain margin ∆. When gain margin is high, the ripple magnitude becomes infinitely large, consequently, the repeater receives multiple peaks of power depending on the number of channels in their frequencies.

Nevertheless, a positive gain margin without AGC or RFE cancellation system cannot exist in practice, because electronics circuits may be damaged before reaching zero gain margin.

In order to characterize this ripple phenomenon, it is also interesting to determine the delay and the attenuation. P (ω) is a periodic function whose period ρ(τ ) can be expressed, from eq. ( 12) such as

ρ(τ ) = 2π τ , (17)
while, attenuation of the channel can be written so that

A(R) = 20 log 10 R/20 -1 10 R/20 + 1 -G. ( 18
)
From eq. ( 17) and eq. ( 16), dualistically, the main parameters of the echo, i.e. ∆ and τ , are given by

         ∆(R) = 20 log 10 R/20 + 1 10 R/20 -1 , τ (ρ) = 2π ρ . (19)
Therefore, as ripple magnitude R and the period ρ can be measured by using frequency sweeping technique, it is possible to measure indirectly ∆ and τ . Therefore it is possible to characterize RFE, which is presented in the next section.

III. SIMULATION

In this part, the theoretical results are compared to simulated ones by using Simulink, where a signal generator block provides a sinusoidal signal with AWGN and both signal and noise parameters are independently adjustable. The echo creator is composed of a repeater (gain, delay, noise) and a channel (attenuation, delay) model in closed-loop configuration. Thus, the first echo is added to the principal signal after τ delay ; the l th echo is added after lτ delay as expressed in eq. ( 4). Echoes are reduced by a A dB factor for every passing through the loop. Hence, the first echo is reduced by A dB and the l th echo is reduced by lA dB. A peak detector block is realized by using the Fast Fourier Transform (FFT) in order to calculate the Power Spectral Density (PSD). The peak detector block provides the value of the maximum (in dB) for any frequency value, when sweeping. Finally, the result exporter block exports simulation parameters and results to a database. Unless otherwise specified, the sampling frequency is set to 100 MHz. The frequency sweep ranges from 20 MHz to 30 MHz with a frequency sweep step of 1 KHz and is realized with a cosine signal without any modulation. The channel attenuation is set to -50 dB. Finally, the other variables are depending on each case. A Gui window is used to control system parameters automatically, sweeping frequency, gain, etc. Fig. 4 shows an example of sweeping frequency, with a frequency range of 200 KHz. It provides a comparison of the maximum of PSD value calculated in theory, in equation (12), and simulation results. A systematic study has been performed for different values of ∆ and τ , showing again a very good match between simulation results and theoretical predictions. It has also been observed that the measures of maximum of PSD are only slightly sensitive to noise, as illustrated in Fig. 4. Concerning the ripple magnitude R, a systematic study in function of the gain margin is presented in Fig. 5 where theoretical results and simulation results for different values of AWGN can be compared. As it can been observed, for a high SNR, errors between simulation results and theory (equation 16) are negligeable and even for low SNR (for instance 30 dB), the errors remain very small. In equation ( 19), it is shown that the gain margin ∆ is depending only on the ripple magnitude R and the loop delay τ is depending only on the period of ripple ρ. By using an appropriate system of data acquisition, as digital acquisition using a Field-Programmable Gate Array (FPGA), R and τ can be estimated using a nonlinear estimating technique as leastsquares [START_REF] Marquardt | An algorithm for least-squares estimation of nonlinear parameters[END_REF], Gauss-Newton or Levenberg-Marquardt [START_REF] Gavin | The Levenberg-Marquardt method for nonlinear least squares curve-fitting problems[END_REF] to cite but a few. Hence ∆ and τ can be calculated using eq. ( 19) with the estimated value of R and ρ.

In the following, the effect of AWGN, the gain margin and the band width on the Relative Error Estimation (REE) are discussed.

Since it has been demonstrated that ripple magnitude are depending on the gain margin value (eq. 16), estimation precision is naturally depending on the gain margin value which is corrupted by the AWGN. Top part of Fig. 6 shows the REE of the estimation of ∆ depending on the gain margin value for different values of SNR. According to the figure, gain margin is best estimated for SNR values more than 30 dB and gain margin values between -10 and -5 dB. Loop delay estimation is optimal for any SNR value according to the lower part of Fig. 6, as the corresponding REE remains lower than 0.14 %, even if the SNR reaches 30 dB.

Sweeping frequency can be a slow operation (with a large frequency range and a small sweeping step). In practice, emitting signals must be in accordance with each country's rules and regulations, this is why the bandwidth effect is studied in this section. Indeed, it is important to reduce to a minimize the range of sweeping. Top part of Fig. 7 shows that REE, when estimating gain margin, is less than 5% for SN R > 30 dB and bandwidth larger than 2 MHz, while bottom part of Fig. 7 shows that REE, when estimating loop delay, is less than 0.5% for SN R > 30 dB and bandwidth larger than 2 MHz. In all cases, SNR values, which corresponds to the quality of the signal, implies a limitation for an accurate estimation of the parameters characterizing RFE. Simulations have also shown that eq. ( 19) is satisfied, i.e, echo parameters are well measured.

CONCLUSION

This paper discussed a method for radio frequency echo detection and measurement by sweeping carrier frequency and by analyzing ripples of the power spectrum density. Mathematical expressions have been established. These expressions have been confirmed by simulation results.

This solution can be used in addition to an echo cancellation system [START_REF] Lee | A radio repeater interference cancellation model for mobile communication systems[END_REF] to improve the efficiency and shut down the echo cancellation system during the absence of radio frequency echo in order to economize energy.

In a simplified model of radio frequency echo, the main parameters are the loop delay (τ ) and the gain margin (∆). These parameters are calculated in theory and confirmed by simulation. Nonlinear estimation methods have been used in this work. It has also been shown the influence of noise (AWGN) and the range of frequency sweeping on the echo parameter estimation.

This solution could be used to detect the presence of echoes and to trigger alarm after leaving specified value in repeaters, so that it could prevent the signal corrupting and system damage against sudden arrivals of radio frequency echo dues to reflections, climate changes, etc.

Indeed, it has been shown, in simulation part, that the RFE parameters are correctly estimated (less than 5%) in noisy environment SNR≈ 30dB and a range less than 2 MHz with 1 KHz of sweeping step. Under these conditions, radio frequency echo parameters could be used to cancel the echo by creating a feedback and learning loop. It would be possible to enhance these results by increasing system performance as noise reduction in repeater, frequency sampling in digital part, estimating method and so on. As a perspective of future work, this study could be extended to the case of multipath or modulated signals.

Fig. 1 .

 1 Fig. 1. a) Spectrum of wide band signal (10 MHz centered on 35 MHz) and narrow band (200 KHz centered on 45 MHz). b) Echo is given for gain margin ∆ = 5 dB and loop delay τ = 0.5 µs.

  Power : P = 10 log(p) Gain : G = 20 log(g) Attenuation : A = 20 log(a) Gain margin : ∆ = 20 log(δ) = G + A.

Fig. 3 .

 3 Fig. 3. Representation of power with strong (top figure) and weak (bottom figure) echo for G = 50 dB and τ = 0.5 a.u. Top figure: ∆ = [-5, -4, -3, -2, -1] dB. Bottom figure: ∆ = [-10, -15, -20, -25, -30] dB.

Fig. 4 .

 4 Fig. 4. The maximum of PSD, with A = -50 dB, G = 45 dB (∆ = 5 dB), τ = 0.5 µs, without noise (dashed curve) and with SN R = 30 dB.

Fig. 5 .

 5 Fig. 5. The ripple magnitude R depending on the gain margin ∆, comparing theory and simulation for different values of SNR.

Fig. 6 .

 6 Fig. 6. Gain margin ∆ (top figure) and Loop delay τ (bottom figure) estimation. Input parameters : ∆ = 10 dB for τ estimation, τ = 5 µs for ∆ estimation.

Fig. 7 .

 7 Fig. 7. Gain margin (∆ = 10 dB) and Loop delay (τ = 5µs) estimation for different values of SNR.

ACKNOWLEDGMENT

The authors would like to acknowledge support from the Conseil Régional de Bourgogne (Regional Council of Burgundy) in the framework of the PARI program and ANRT agency (French National Agency for Research and Development) in the framework of CIFRE program.


Santosh Tirunagari
email: santosh.tirunagari@surrey.ac.uk

Samaneh Kouchaki
email: samaneh.kouchaki@manchester.ac.uk

Norman Poh
email: n.poh@surrey.ac.uk

Miroslaw Bober
email: m.bober@surrey.ac.uk

David Windridge
email: d.windridge@mdx.ac.uk.

Dynamic Mode Decomposition for Univariate Time Series: Analysing Trends and Forecasting

Keywords: Time series modelling, dynamic mode decomposition, DMD, forecast, prediction, trends, noise

This study introduces the method of Dynamic Mode Decomposition (DMD) for analysing univariate time series by forecasting as well as extracting trends and frequencies. The key advantage of DMD is its data-driven nature which does not rely on any prior assumptions (like Singular Spectrum Analysis (SSA)) except the inherent dynamics which are captured over time. Indeed, this study will show that the DMD eigenvalues with frequencies that are closer to the origin in the complex plane capture the trends in the time series. Moreover, the temporal evolution of the DMD modes, which is preserved via the Vandermonde matrix, can be used to reconstruct the desired components and perform forecasting at the same time. The results at various noise levels on simulated data suggests that DMD is a promising approach to modelling a time series with a noisy structure. Although these properties are not new in the DMD literature, the novel contributions of this paper are in making the method of DMD work for univariate time series through a four staged pipeline. Thus, this is the first work that shows DMD can be used for modelling, predicting, and forecasting a univariate time series.

I. INTRODUCTION

Time series analysis includes extracting meaningful information like filtering trends, and forecasting future values [START_REF] Golyandina | Singular Spectrum Analysis for time series[END_REF], [START_REF] Box | Time series analysis: forecasting and control[END_REF]. It has applications in many areas, such as quantitative finance and biomedical engineering [START_REF] Sanei | Singular spectrum analysis of biomedical signals[END_REF], [START_REF] Mohammadi | Improving time-frequency domain sleep {EEG} classification via singular spectrum analysis[END_REF], [START_REF] Kouchaki | Tensor based singular spectrum analysis for automatic scoring of sleep eeg[END_REF]. This signal processing problem is often tackled by decomposing the data into several components in the frequency domain, thus decoupling signals or trends from noise. To this end, matrix decomposition techniques such as singular spectrum analysis (SSA) [START_REF] Sanei | Singular spectrum analysis of biomedical signals[END_REF], [START_REF] Elsner | Singular spectrum analysis: a new tool in time series analysis[END_REF] is often used. Analogous to SSA, this paper presents a data-driven approach, Dynamic Mode Decomposition (DMD).

The method of DMD was introduced originally in the area of Computational Fluid Dynamics (CFD) for extracting coherent structures from spatio-temporal complex fluid flow data [START_REF] Schmid | Dynamic mode decomposition of numerical and experimental data[END_REF]. DMD takes in time series data and computes a set of modes, each of which is associated with a complex eigenvalue. DMD analysis is closely associated with spectral analysis of the Koopman operator, which provides linear but infinite-dimensional representation of nonlinear dynamical systems [START_REF] Tu | On dynamic mode decomposition: theory and applications[END_REF]. Therefore, by using DMD a nonlinear system could be described by a superposition of modes whose dynamics are governed by the eigenvalues. The key advantage of DMD is its data-driven nature which does not rely on any prior assumptions except the inherent dynamics which are observed over time. Its capability for extracting relevant modes from complex fluid flows [START_REF] Schmid | Dynamic mode decomposition of numerical and experimental data[END_REF], [START_REF] Schmid | Dynamic mode decomposition and proper orthogonal decomposition of flow in a lid-driven cylindrical cavity[END_REF], [START_REF] Schmid | Applications of the dynamic mode decomposition[END_REF], [START_REF] Tirunagari | Analysis of proper orthogonal decomposition and dynamic mode decomposition on les of subsonic jets[END_REF] has seen significant application across multiple fields, including computer vision, for robustly separating video frames into a background model and multiple foreground objects [START_REF] Grosek | Dynamic mode decomposition for real-time background/foreground separation in video[END_REF], [START_REF] Tirunagari | International Conference on Image, Vision and Computing (ICIVC)[END_REF], [START_REF] Erichson | Randomized low-rank dynamic mode decomposition for motion detection[END_REF]. In biometrics for detecting spoofed samples [START_REF] Tirunagari | Detection of face spoofing using visual dynamics[END_REF], [START_REF] Tirunagari | Windowed dmd as a microtexture descriptor for finger vein counter-spoofing in biometrics[END_REF] from the live ones. In robotics and neuroscience, DMD has been used to estimate perturbation in human robot interactions [START_REF] Berger | Dynamic mode decomposition for perturbation estimation in human robot interaction[END_REF] and to extract coherent patterns in large-scale neural recordings [START_REF] Brunton | Extracting spatial-temporal coherent patterns in large-scale neural recordings using dynamic mode decomposition[END_REF].

All the aforementioned DMD applications consider a data matrix as an input. For instance, to analyse an image sequence (a video) using DMD, the input data is converted into column vectors, each of which representing a single frame of image. However, being a fundamentally multidimensional method, DMD has never been applied to a univariate time series. Therefore, our contribution is to introduce the method of DMD through a four staged approach for decomposing a univariate time series into a number of interpretable elements in terms of noise, trend and harmonics. Moreover, time series forecasting using DMD is also presented in this study.

The methodological pipeline consists of four stages similar to that of SSA. However, it is different from SSA in the formulation. The beauty of DMD is that it captures frequency components that describe the growth/decay and oscillations in the data whilst at the same time it can forecast time series during the reconstruction stage.

In the remainder of this paper, the DMD methodology is described in Section II and SSA is described in Section III. Section IV presents the datasets used in this study. Experiments and results are presented in Section V. Finally, conclusions with discussions are drawn in Section VI.

II. METHODOLOGY

This section formulates DMD for univariate time series analysis. Our algorithm consists of four stages: (1) embedding, (2) dynamic mode decomposition, (3) reconstruction and forecasting, and (4) diagonal averaging. Each step is explained in more detail in the following subsections:

A. Embedding

Embedding can be considered as a mapping that converts a univariate time series data

x = [x 1 , x 2 , • • • , x T ] of length T into a multidimensional series X = [x 1 , x2 , ..., xL ], of size K × L. The vector xn = [x n , x n+1 , • • • , x n+K-1 ] T ∈ R K is called a K-lagged vector, and L ∈ [2 ≤ L ≤ T 2 ]
is the length of window whereas K = T -L + 1, is the number of segments. Thus, the time series is decomposed into K overlapping segments of length L. This produces a trajectory matrix X ∈ R K×L . Unfolding the matrix X with elements K × L, we have:

X =     x 1 x 2 ... x L x 2 x 3 ... x L+1 . . . . . . . . . . . . x K x K+1 ... x K+L-1     . (1)
Note that X is a Hankel matrix, because its anti-diagonal elements (m + n) = const, are equal. Where, m,n are the row and column indices of the Hankel matrix1

B. Dynamic Mode Decomposition

The K-lagged vectors xn ∈ R K are evenly spaced by δt = K. Therefore, it can be justified that a mapping A exists between consecutive lagged vectors. Formally, we can say that these vectors span the krylov subspace [START_REF] Krylov | On the numerical solution of the equation by which in technical questions frequencies of small oscillations of material systems are determined[END_REF], [START_REF] Saad | Krylov subspace methods for solving large unsymmetric linear systems[END_REF], [START_REF] Ruhe | Rational krylov sequence methods for eigenvalue computation[END_REF], i.e.,

P = [x 1 , Ax 1 , A 2 x1 , A 3 x1 , • • • , A N -1 x1 ] (2)
for an unknown matrix A, this formulation (Eq. ( 2)), essentially reconstructs the data X at any instance

l ∈ {1, • • • , L - 1, • • • L -1 + F } after the initial vector x1 .
Here L is the number of columns in the trajectory matrix X and F is the number of future points. This will be discussed in Section II-C. Eq. ( 2) can be rewritten by introducing two data matrices

P 2 ≡ [x 2 , x3 , • • • , xL ] and P 1 ≡ [x 1 , x2 , • • • , xL-1 ]
in the following way:

P 2 = AP 1 . (3)
The mapping matrix (also known as Koopman operator) A is responsible for capturing the dynamics within these vectors. The sizes of the matrices P 2 and P 1 are K × (L -1) each. Therefore, the size of the unknown matrix A would be K × K. For (K L), unfortunately, solving for A will be computationally very expensive. Moreover, our problem of interest is to find the temporal evolution which are essentially captured in terms of column vectors. Therefore, it is necessary to represent these column vectors as a linear combination. Our assumption that these vectors form a Krylov span, allows us to introduce H (as a starting point of the standard Arnoldi iteration),

AP 1 ≈ P 1 H. (4)
Here, H is a companion matrix, also known as a shifting matrix, that simply shifts vectors 1 through L -2. In essence, H approximates the last vector L -1 by linearly combining the previous L -2 vectors, i.e., xL-1 = c 0 x1 + ...

+ c L-2 x1 .
H requires the storage of (L -1) × (L -1) data matrix and is significantly fewer in dimension than A.

H =       0 0 . . . 0 c 0 1 0 . . . 0 c 1 0 . . . . . . . . . . . . 0 0 1 0 c L-2 0 0 0 1 c L-1       . (5)
Thus, for the last K lagged vector xL-1 , where L is significantly fewer in dimensionality of A (K × K), one can write P 2 as a linear combination of the previous vectors. Consistent with Eqs. ( 3) and (4), we then have:

P 2 = P 1 H + r.e * L-1 , (6)
where, r is the residual error, (.) * is the conjugate transpose operator, e L-1 is the (L -1) unit vector. From Eqs. ( 4) and ( 6), we have AP 1 ≈ P 2 ≈ P 1 H. Although we can solve H using the LU [START_REF] Tirunagari | Detection of face spoofing using visual dynamics[END_REF] or the QR decomposition [START_REF] Schmid | Applications of the dynamic mode decomposition[END_REF], [START_REF] Lawson | Solving least squares problems[END_REF], a more robust solution can be achieved using SVD [START_REF] Schmid | Applications of the dynamic mode decomposition[END_REF] on P 1 in Eq. ( 3), to obtain U, Σ and V * matrices.

P 2 = AUΣV * . ( 7
)
∵ P 1 = d i=1 Σ i u i v * i ,
here, Σ i is the i th singular value of P 1 , u i and v * i are the corresponding left and right singular vectors respectively; and d is the total number of singular values. Rearranging Eq. ( 7), we obtain the the full-rank matrix A,

A = P 2 VΣ -1 U * . ( 8
)
Since the eigenvalue analysis is agnostic to any linear projection; so solving the eigen problem of H is easier than that of solving for A directly. Moreover, the associated eigenvectors of H provide the coefficients for the linear combination that is necessary to express the dynamics within the time series basis.

Hω = σω, (9)
where, ω are the eigenvectors and σ a diagonal matrix containing the corresponding eigenvalues of H matrix. The eigenvalues of H approximate some of the eigenvalues of the full system A [13], we then have:

AU ≈ U H, AU ≈ Uωσω -1 , A(Uω) ≈ (Uω)σ. ( 10
)
Therefore H is determined on the subspace spanned by the orthogonal singular basis vectors U obtained via P 1 ,

H = U * (A)U, H = U * (P 2 VΣ -1 U * )U, (11)
which can be rewritten as:

H = U * P 2 VΣ -1 . ( 12
)
Here U * ∈ C (L-1)×K and V ∈ C (L-1)×(L-1) are the conjugate transpose of U and V * , respectively; and Σ -1 ∈ C (L-1)×(L-1) denotes the inverse of the singular values Σ.

By replacing Ψ = Uω in Eq. ( 10) i.e., A(Ψ) ≈ (Ψ)σ, we obtain the dynamic modes Ψ. ∵ U = P 2 VΣ -1 , therefore, we have:

Ψ = P 2 VΣ -1 ω (13)
The complex eigenvalues σ contain growth/decay rates and frequencies of the corresponding DMD modes [START_REF] Schmid | Dynamic mode decomposition and proper orthogonal decomposition of flow in a lid-driven cylindrical cavity[END_REF], [START_REF] Schmid | Dynamic mode decomposition of numerical and experimental data[END_REF]. If σ j are the diagonal elements of σ from Eq. ( 9), the temporal behaviour of the DMD modes is then formed via Vandermonde matrix V and , which raises its column vector to the appropriate power. V and (f ) with (L-1)×(f +1) elements will be defined by the following:

V and (f ) =      1 σ 1 1 σ 2 1 ... σ f 1 1 σ 1 2 σ 2 2 ... σ f 2 . . . . . . . . . . . . . . . 1 σ 1 L-1 σ 2 L-1 ... σ f L-1      , (14)
V and (L) is a standard Vandermonde matrix for reconstruction but if f > L, this is used for forecasting, which will be discussed in section II-C.

DMD modes with frequencies µ j is defined by:

µ j = ln(σ j ) δt , ( 15
)
where δt is the lag between the vectors and may be considered to consist of K time steps. The real part of µ j regulates the growth or decay of the DMD modes, while the imaginary part of µ j drives oscillations in the DMD modes. The frequencies near the origin can be interpreted as the trend of the time series, and the frequencies bounded away from the origin contain harmonics of the time series. Specifically, the trend in the time series has an associated frequency at the origin of the complex plane with µ j ≈ 0.

C. Reconstruction and Forecasting

The DMD reconstruction of the data X at any instance f after the initial vector x1 is given by,

X = Ψ(V and (f ) • b), ∀ f ∈ {1, 2, • • • , F }. ( 16
)
Where V and (f ) = σ f -1 and F is the number of future points that are to be forecast. The vector b ≈ Ψ -1 x 1 contains the initial amplitudes for the dynamic modes. Operator (•) is the point-wise (element-wise) multiplication of b with every column of V and (f ). The first lagged vector x1 of the data X reduces to

x1 = Ψb, (17)
and b can be obtained through the Moore-Penrose pseudo

inverse, b = (Ψ * Ψ) -1 Ψ * x1 . ( 18
)
The significance of the DMD modes Ψ are determined by their scaling factors b when it comes to approximating the entire data sequence P 2 during reconstruction and forecasting.

Since the slow varying frequencies µ j capture trends in the time series, it is desirable to reconstruct the trajectory matrix with selected DMD triples. A DMD triple consists of a DMD eigenvalue, its corresponding DMD mode and an associated scaling factor, i.e., amplitude.

X = Ψ {S} (V and{S} (f ) • b {S} ), S ⊂ {1, ..., L -1}. (19)
V and{S} (f ) is a subset of V and (f ) formed by selecting only the appropriate rows which correspond to the selected DMD eigenvalues. V and{S} (f ) will have the dimension of S × (f + 1) and b {S} has the dimension of (S × 1) for the selected S eigenvalues.

D. Hankelisation and Diagonal Averaging

Finally, the reconstructed trajectory matrix X is Hankelised, in order to ensure that the anti-diagonal elements are equalised. The Hankelisation operator H for an K×F matrix X is defined as

H X = X =     x 1 x 2 ... x F x 2 x 3 ... x F +1 . . . . . . . . . . . . x K x K+1 ... x K+F -1     . (20)
The desired or reconstructed time-series is given by

( x k ) = 1/num(D k ) m,n∈D k Xm,n .
Here k is the index of the time series, (.) is the real part of x k , num(.) is the number of combinations of (m, n), m

+ n = k + 1 and D k is given by {(m, n) : 1 ≤ m ≤ K, 1 ≤ n ≤ F, m+n = k+1}.
The methodology, as summarised via utilising Eq. ( 19), can be used to decompose data into a number of desired subcomponents including trend, noise and future components, but at the expense of introducing a time stamp delay. This is because, P 2 has been used for projection from Eq. ( 13), which does not contain information about the first data point. Therefore, this study recommends that the first column vector of the trajectory should be repeated.

III. SINGULAR SPECTRUM ANALYSIS

SSA is another powerful method for analysing real-valued time series [START_REF] Golyandina | Analysis of time series structure: SSA and related techniques[END_REF]. The initial idea of SSA is associated with the work of Broomhead [START_REF] Broomhead | Extracting qualitative dynamics from experimental data[END_REF], [START_REF] Broomhead | On the qualitative analysis of experimental dynamical systems[END_REF]. It combines the multivariate statistics, classical time series analysis, dynamical systems, and signal processing.

SSA is becoming an effective method in various areas such as economics [START_REF] Golyandina | Analysis of time series structure: SSA and related techniques[END_REF], [START_REF] Alexandrov | The automatic extraction of time series trend and periodical components with the help of the Caterpillar-SSA approach[END_REF] and biomedical engineering [START_REF] Sanei | A new adaptive line enhancer based on singular spectrum analysis[END_REF], [START_REF] Ghaderi | Localizing heart sounds in respiratory signals using singular spectrum analysis[END_REF], [START_REF] Sanei | An adaptive singular spectrum analysis approach to murmur detection from heart sounds[END_REF]. Basically, SSA decomposes a data into a number of interpretable elements with different subspaces, such as noise and trend [START_REF] Kouchaki | Supervised single channel source separation of eeg signals[END_REF], and can be used for any time series with complex structure [START_REF] Alexandrov | The automatic extraction of time series trend and periodical components with the help of the Caterpillar-SSA approach[END_REF]. For instance, for decades, SSA has been used for both trend detection and prediction in financial data [START_REF] Golyandina | Analysis of time series structure: SSA and related techniques[END_REF], [START_REF] Alexandrov | The automatic extraction of time series trend and periodical components with the help of the Caterpillar-SSA approach[END_REF]. SSA can be used in several applications including trends fitting, extraction of cycles with various periods and amplitudes, smoothing, and finding some structures in short time series.

Ordinary SSA is a subspace decomposition algorithm with four stages similar to that of the proposed method but with different formulation in the second and third stages. SSA decomposes the resulted matrix of embedding step using SVD:

C = d i=1 C i = d i=1 √ σ i u i v i T , (21)
where σ i is the ith eigenvalue of covariance matrix C = XX T , u i is the corresponding eigenvector, d is the total number of eigenvalues, and

v i = C T u i / √ σ i .
In the reconstruction step, first the elementary matrices from the previous stage are grouped into several sub-matrices

C = Q q=1 Ĉq , (22)
where Q determines the total number of groups, index q refers to qth subgroup of eigenvalues, and Xq indicates the sum of X i within group q. Finally, the desired subgroup is selected and transformed into a Hankel matrix. Later, the average of the anti diagonal elements form the reconstructed desired time series.

SSA has recurrent and vector based forecasting techniques:

A. SSA Recurrent Forecasting

If the desired group of eigentriples is considered as I and the corresponding ith eigenvector as u Ii , then

r = 1 1 -υ 2 i∈I π i u Ii ( 23
)
where r = (a l1-1 , . . . , a 1 ) T , u Ii and π i respectively show all except for the last and the last column of u Ii , and υ2 = i∈I π 2 i . Consequently, the recurrent forecasting can be formulated as

x i =    x i i = 1, . . . , T l1-1 j=1 a j x i-j i = T + 1, . . . , T + F ( 24
)
where F is the number of forecasting points. The reconstructed terms and Eq. ( 23) are used to forecast the value of new points. The linear operator ϕ can be defined as

ϕx = x r T x , (25)
where x shows K -1 components of vector x. Therefore, trajectory matrix of the time series x of length T + F can be written as the vector representation of Eq. ( 24)

z i = x i i = 1, 1, . . . , L ϕz i-1 i = L + 1, . . . , T + F . (26)
IV. DATASET

A. Real Data

To demonstrate and validate the methodology, publicly available Box & Jenkins airline data [START_REF] Box | Time series analysis: forecasting and control[END_REF] have been used. This consists of the monthly total of international airline passengers from the year 1949 to 1960. For comparing with SSA, BJ sales and UK gas data are also considered. All of the aforementioned datasets have been obtained from Github repository for R time series datasets 2 . These datasets then have been transformed into the logscale (y-axis) and the time (x-axis) is considered in increments of 1.

B. Simulated Data

The simulated data are a mixture of two signals

x(t) = as(t) + βe(t), (27)
where s(t) shows the desired signal and a is the mixing vector.

x(t) is the input or mixture being examined, e(t) indicates the unwanted signal or noise, and β is the noise level.

Signal-to-noise ratio (SNR) in terms of mean square (MS) is used as a measure for noise level which is adjusted by changing

β SNR = MS(s(t))/β 2 . ( 28
)
To evaluate the simulation performance, the evaluation criteria in terms of mean square error (MSE) is considered as

MSE = MS(s(t) -ŝ(t)) MS(s(t)) , ( 29
)
where ŝ(t) is the estimated source.

V. EXPERIMENTS AND RESULTS

The real time series containing airline data is transformed to a trajectory matrix with K-lagged vectors. The parameter that needs to be predetermined is the window length L because our conjecture is that an appropriate choice of this can produce a better result.

A. Effect of Window Length L

Fig. 1 (a) shows our experimental results with window lengths L for reconstructing the original time series. Our result shows that the reconstruction error for L ∈ {12, ..., 66} gave a stable performance with relatively low reconstruction error. If L is relatively small, it can cause the trajectory matrix X to miss the trends in the data especially when dealing with seasonal components. Furthermore, if L is closer to T 2 , where T is the length of the time series, the decomposition of the data becomes more detailed. Therefore, the selection of L can neither be close to 2 nor T 2 . Therefore, this study, recommends choosing L between 8% to 45% of the length of the time series.

In our experiments we arbitrarily considered L to be 25% of the length of the time series, i.e., L = 35.

B. DMD for Extracting Trends

The logarithmic values of DMD eigenvalues give us the frequencies, the real part captures the growth/decay and imaginary part, the oscillations within the time series. The frequencies that are closer to origin are the slow varying components (for example, eigenvalues {18, 19, 20, 21, 22} in Fig. 1 (b) which are vertically closer to the X-axis) that essentially capture the trends. The ones which are further away from the origin capture various frequency growth/decay rates. In order to establish the order in the DMD, the absolute frequencies of the DMD eigenvalues are sorted in the ascending order. The absolute value of the frequencies associated with each of the eigenvalues are shown in Fig. 3

(b).

Reconstruction of the original time series with the aforementioned slow varying DMD triples, i.e., with {18, 19, 20} in Fig. 1 (b) (or the first three DMD triples) gave us the trend (trend-1) as shown in Fig. 2 (a).

Reconstruction with other DMD triples that are closer to the origin, i.e., {18, 19, 20, • • • , 26} gave a smoother version (trend-2) of the time series (Fig. 2 (a)).

(Ψj, σj) As we reconstruct the data with increasing DMD triples, the reconstruction mean square error (MSE) decreases (Fig. 3 (a)). A complete reconstruction of the original time series with all of the DMD triples produces the original time series as shown in Fig. 2 (b) with a reconstruction MSE of 0.0006.

C. DMD for Time Series Forecasting

To demonstrate the capability of DMD for forecasting, 124 out of the total 144 datapoints are considered as the training set, leaving the remaining 20 points as the test set. The window length L considered comprises of 25% of the data i.e., L = 31. The forecast trends in the data are shown in Fig. 4 (a). The MSE recorded for the forecast trend-1 was 0.0007 and for the smoother trend-2 was 0.0080 when evaluated against the reconstructed trend. The result of the forecast is shown in Fig. 4 (b) with a MSE of 0.0090 when compared to the original data.

D. Comparison to SSA 1) On Simulated Data with Effect of Noise: In order to compare with SSA, the performance of both the methods should be investigated in the presence of noise. For this purpose, simulated signal consisting of two sinusoids with different frequencies and powers are considered, which were then mixed with a Gaussian noise at various signal to noise ratio (SNR) levels {-10, -9, • • • , 19, 20}dB.

2) Selection of the triples: The proposed method using the first 9 DMD triples and SSA with first 4 singular triples were applied as they those parameters gave the least reconstruction MSE across all the noise levels on this simulated data as shown in Figure5. 3) Robustness to Noise: It can be noticed from Fig. 6 that, as the SNR increases, the MSE in reconstruction (with respect to the clean signal) decreases in both of the methods. The average time taken for DMD across all the noise levels was 3.76, while for SSA was 46.94 seconds 3 . This result shows that DMD is computationally more efficient than SSA and is consequently a promising approach for analysing a univariate time series in the presence of noise. Fig. 7 shows the original simulated signal and the noisy signal at SNR = -2dB. The reconstruction along with forecasting of 750 points, using the first 5 DMD triples and 4 singular triples is shown. The MSE recorded after reconstruction and forecast was 0.0205 and for SSA was 0.0241. The time taken for DMD and SSA was 0.17 and 23.23 seconds respectively. I. These results show that DMD based univariate time series analysis is superior to SSA particularly when forecasting the time series.

SNR (dB)

VI. CONCLUSION AND DISCUSSIONS

In this study we proposed the method of DMD in a fourstage pipeline to decompose a univariate time series into a number of interpretable elements with different components, such as noise and underlying data trend. In addition, this fourstage DMD pipeline was also shown to be used for forecasting a time series. DMD eigenvalues with frequencies that are closer to the origin in the complex plane capture the trends in the time series. Moreover, the temporal evolution of the DMD modes, which is preserved via the Vandermonde matrix, can be used to reconstruct the desired components and perform forecasting at the same time. Window length L needs to be predetermined and we recommend choosing L between 8% to 45% of the length of the time series. In reconstruction, to compensate one time stamp delay, we recommend that the first column vector of the trajectory matrix be repeated once. The experimental results at various noise levels on simulated data suggests that DMD is a promising approach to modelling a time series with a noisy structure. DMD was compared to SSA, a popular method employed extensively in biomedical signal processing applications for separating artefacts and extracting features from electromyography (EMG), electrocardiogram (ECG) and electroencephalogram (EEG) signals [START_REF] Sanei | A new adaptive line enhancer based on singular spectrum analysis[END_REF].

Thus DMD can potentially be applied to analyse biomedical signals, especially EEG, where brain rhythms manifest themselves as narrow frequency band components. As a good example, sleep is a dynamic process which consists of different stages with different neural activity levels. Each stage is characterised by a distinct set of physiological and neurological features and dominant frequency band. SSA, however, does not exploit the narrowband property of these kind of cyclic data since it requires manual grouping or optimisation [START_REF] Sanei | An adaptive singular spectrum analysis approach to murmur detection from heart sounds[END_REF] to detect the desired features. DMD on the other hand due to its intrinsic nature of complex eigenvalues that capture frequency bands, can help find cycling components automatically. Therefore, in such studies, DMD might capture clues for chaotic subcomponents/trends and noise. These subcomponents later can be used as feature vectors for e.g. in sleep stage classification, epilepsy detection, and seizure predictions (for e.g., using SSA in [START_REF] Sanei | A new adaptive line enhancer based on singular spectrum analysis[END_REF], [START_REF] Sanei | Singular spectrum analysis of biomedical signals[END_REF], [START_REF] Mohammadi | Improving time-frequency domain sleep {EEG} classification via singular spectrum analysis[END_REF], [START_REF] Kouchaki | Tensor based singular spectrum analysis for automatic scoring of sleep eeg[END_REF]. Moreover, our experiments showed that DMD is computationally more efficient than SSA as well as robust in the presence of added white Gaussian noise.

Fig. 1 .

 1 Fig. 1. (a) Effect of window length L over the reconstruction performance. (b) DMD eigenvalue plot, showing eigenvalues {18, 19, 20, 21, 22} that are vertically closer to the X-axis

Fig. 2 .

 2 Fig. 2. (a) Trend analysis of the timeseries, revealing two trends (i) trend1: linear trend and (ii) periodic trend. (b) A complete reconstruction.

Fig. 3 .

 3 Fig. 3. (a) DMD triples vs reconstruction error. (b) Corresponding absolute frequencies.

Fig. 4 .

 4 Fig. 4. (a) Trend forecast. (b) A complete forecast of the timeseries.

Fig. 5 .

 5 Fig. 5. MSE while reconstruction of the original time series in the presence of noise at -2dB, using first t (a) DMD triples, (b) SVD triples.

Fig. 6 .Fig. 7 .

 67 Fig.6. Effect of noise in reconstruction at various SNRs across 500 iterations.

Fig. 8 .

 8 Fig. 8. Comparison of DMD and SSA on (a) Airline passengers data, (b) BJ sales data) and (c) UK gas data.

TABLE I .

 I THE PARAMETERS THAT HAVE BEEN USED FOR ANALYSING THE REAL DATA USING DMD AND SSA. HERE PARAMETERS, T AND L REPRESENT THE TOTAL LENGTH OF THE TIME SERIES AND LENGTH OF THE WINDOW CONSIDERED FOR THE ANALYSIS RESPECTIVELY. 'REC' DENOTES RECONSTRUCTION AND 'FOR' DENOTES FORECASTING. MSE DENOTES THE MEAN SQUARE ERROR BETWEEN THE RECONSTRUCTED/FORECAST TIME SERIES AND THE ORIGINAL TIMESERIES.

	Dataset	Parameters			DMD Triples			SSA Triples			MSE in Rec ,For
		T		L		Trend		Period		Trend		Period		DMD	SSA
		Rec	For	Rec	For	Rec	For	Rec	For	Rec	For	Rec	For
	Airline Passengers	144	130	36	33	2	2	7	9	1	2	3	3	0.0010	0.0004
	BJSales	150	135	38	34	2	2	7	7	1	1	7	7	0.00006	0.0001
	UKgas	97	87	24	22	2	2	12	12	1	2	2	2	0.0088	0.0142

The authors in[START_REF] Brunton | Extracting spatial-temporal coherent patterns in large-scale neural recordings using dynamic mode decomposition[END_REF] use block-Hankelisation approach on the multidimensional data. Our proposed algorithm is especially for univariate data.

https://github.com/vincentarelbundock/Rdatasets

On a Dell PowerEdge R715 AMD -2-CPU/24-Core -128GB RAM Running Precise (Ubuntu 12.04.5).

ACKNOWLEDGEMENT

The funding for this work has been provided by the Department of Computer Science and the Centre for Vision, Speech and Signal Processing (CVSSP) -University of Surrey. 'S.T' and 'N.P' have benefited from the Medical Research Council (MRC) funded project "Modelling the Progression of Chronic Kidney Disease" under the grant number R/M023281/1. The details of the project are available at www.modellingCKD.org.