nonlinearity/spectrum/impl_level_vs_time.py

+"""
+Created on: 11 Jun 2025 - 13:56
+Author: jan.reimes@head-acoustics.com
+"""
+import numpy as np
+from numpy.typing import NDArray
+from functools import partial
+
+from .misc import get_value_representation
+from .types import ISpectrumVsTime, ILevelVsTime,Weighting, ValueRepresentation, Ax
+
+
+class LevelVsTime(ILevelVsTime):
+
+    def __init__(self, spec_vs_time: ISpectrumVsTime):
+        self._spec_vs_time = spec_vs_time
+
+    @property
+    def time(self) -> NDArray:
+        return self._spec_vs_time.time
+
+    @property
+    def fs_out(self) -> float:
+        return 1/(self._spec_vs_time.time[1] - self._spec_vs_time.time[0]) if self._spec_vs_time.time.size > 1 else 0.0
+
+    def get_values(self, representation: ValueRepresentation = 'dB', weighting: Weighting = None) -> NDArray:
+        lvt = self._spec_vs_time.get_values('pow', weighting, value_transform=partial(np.sum, axis=Ax.FREQ),)
+
+        # always include window correction for level vs. time:
+        lvt *= self._spec_vs_time.spectrum_args.get_window_correction('pow')
+
+        return get_value_representation(lvt, representation, 'pow')
+
+    def get_active_frames(self, inactivity_fraction: float, weighting: Weighting = None) -> NDArray[np.bool_]:
+        lxx = self.get_values('dB', weighting)
+        activity_threshold_db = np.percentile(lxx, 100 * inactivity_fraction)
+        return lxx >= activity_threshold_db
+
+
+if __name__ == "__main__":
+    pass
\ No newline at end of file

nonlinearity/spectrum/misc.py

+"""
+Created on: 04 Jun 2025 - 11:44
+Author: jan.reimes@head-acoustics.com
+
+Miscellaneous functions for spectrum calculations. Avoid circular imports!
+"""
+
+from typing import Optional
+import numpy as np
+from numpy.typing import NDArray
+
+from .types import ValueRepresentation
+
+
+def to_db(x: NDArray, is_rms: bool, use_soft_minclip: bool = True, abs_min_value_db: Optional[float] = -100.0) -> NDArray:
+    """
+    Convert a linear scale array to dB scale.
+    Args:
+        x (NDArray): Input array.
+        in_db (bool): If True, convert to dB scale.
+
+    Returns:
+        NDArray: Converted array in dB scale if in_db is True, otherwise unchanged.
+    """
+    # Convert to dB scale
+    a = 20 if is_rms else 10
+    x = np.abs(x)
+
+    if use_soft_minclip:
+        # use smallest valid value in x as minimum
+        if (ii := (x <= 0)).any():
+            if not ii.all():
+                x[ii] = x[~ii].min()
+
+    if abs_min_value_db is not None:
+        # use absolute minimum value
+        x = np.maximum(x, 10 ** (abs_min_value_db / a))
+
+    # final check for negative or zero values - should not happen after applying soft minclip and abs_min_value_db?
+    if (x <= 0).any():
+        raise ValueError("All values in x are negative or zero, cannot convert to dB scale.")
+
+    return a * np.log10(x)
+
+def get_value_representation(x: NDArray, output_representation: ValueRepresentation, input_representation: ValueRepresentation, **to_db_kwargs) -> NDArray:
+    """
+    Get the representation of the input array based on the specified representation type.
+
+    Args:
+        x (NDArray): Input array.
+        rep (ValueRepresentation): Representation type ('rms', 'pow', 'dB').
+
+    Returns:
+        NDArray: The input array transformed to the specified representation.
+    """
+
+    r_in = input_representation.lower()
+    r_out = output_representation.lower()
+
+    for r in (r_in, r_out):
+        if r not in ('rms', 'pow', 'db'):
+            raise ValueError(f"Invalid representation: {r}. Must be one of 'rms', 'pow', or 'dB'.")
+
+    # easy case
+    if r_in == r_out:
+        return x
+
+    # convert if needed
+    if r_out == 'db':
+        return to_db(x, is_rms=r_in == 'rms', **to_db_kwargs)
+    elif r_out == 'rms':
+        if r_in == 'pow':
+            return np.sqrt(x)
+        elif r_in == 'db':
+            return np.sqrt(10 ** (x / 20.0))
+        else:
+            raise ValueError(f"Cannot convert from {r_in} to {r_out}.")
+
+    elif r == 'pow':
+        if r_in == 'rms':
+            return np.square(np.abs(x))
+        elif r_in == 'db':
+            return 10 ** (x / 10.0)
+        else:
+            raise ValueError(f"Cannot convert from {r_in} to {r_out}.")
+    else:
+        raise ValueError(f"Unknown representation: {r}. Must be one of 'rms', 'pow', or 'dB'.")
+
+
+
+if __name__ == "__main__":
+    pass
\ No newline at end of file

nonlinearity/utils/misc.py → nonlinearity/spectrum/spectral_args.py

-# -*- coding: utf-8 -*-
 """
-Created on Sep 27 2024 16:19
+Created on: 06 Jun 2025 - 17:14
+Author: jan.reimes@head-acoustics.com
 
-@author: Jan.Reimes
+This module provides configuration and validation for spectral analysis parameters.
+Contains dataclasses and utilities for handling FFT settings, windowing options,
+and frequency range specifications used across spectrum analysis functions.
 """
-
+from dataclasses import dataclass, asdict
+from typing import Optional, Tuple, Literal, Self, Dict, Any
 import numpy as np
 from numpy.typing import NDArray
-from typing import Union, Iterable, Dict, Any, Optional
 from ensure import Check
 from scipy.signal import get_window
 
+from .types import ValueRepresentation
 
-def interpolate_spectrum(spec: NDArray, freq: NDArray, values_to_interpolate: Union[float, Iterable[float]],
-                         log_frequencies: bool = True, axis: int = 0) -> NDArray:
-
-    # linearly interpolates data points with value(s) given in values_to_interpolate
-    Check(len(spec.shape)).equals(2).or_raise(ValueError, f"Only 2D spectrum supported, but got {spec.shape}")
-
-    # move channel axis to 0
-    if axis != 0:
-        spec = np.moveaxis(spec, axis, 0)
-
-    values_to_interpolate = np.atleast_1d(values_to_interpolate)
 
-    if log_frequencies:
-        freq = np.log2(freq)
-
-    # iterate over all channels
-    for i in range(spec.shape[0]):
-        idx = np.zeros(spec.shape[1], dtype=bool)
-
-        for v in values_to_interpolate:
-            idx = idx | (spec[i] <= v)
+@dataclass(frozen=True)
+class SpectralArguments:
+    """
+    Class to hold spectral arguments for FFT and band spectrum calculations.
+    """
+    fs: int|float
+    window: Optional[str|NDArray] = 'hann'
+    nperseg: int = 8192
+    nfft: Optional[int] = None
+    fmin: float = 20.0
+    fmax: float = 20_000.0
+    overlap: float = 0.75
+    mode: Literal['psd', 'stft'] = 'psd'
+    scaling: Literal['density', 'spectrum'] = 'spectrum'
+    padded: bool = True
+
+    def to_dict(self) -> dict:
+        """
+        Convert the SpectralArguments instance to a dictionary, excluding None values and empty lists/tuples.
+        """
+        return {k: v for k, v in asdict(self).items() if v is not None and not isinstance(v, (list, tuple))}
+
+    @classmethod
+    def from_kwargs(cls, **kwargs: dict) -> Tuple[Self, Dict[str, Any]]:
+        # create a new instance of SpectralArguments from (arbitrary) keyword arguments, returning also the remaining/unused kwargs
+        fields = {}
+        remaining_kwargs = kwargs.copy()
+        for k, v in kwargs.items():
+            if k in cls.__match_args__:
+                fields[k] = v
+                remaining_kwargs.pop(k)
+
+        return cls(**fields), remaining_kwargs
+
+    def __post_init__(self):
+        if not (0 <= self.overlap < 1):
+            raise ValueError("Overlap must be in the range [0, 1).")
+        if self.fs <= 0:
+            raise ValueError("Sampling frequency (fs) must be positive.")
+        if self.nfft is None:
+            # If nfft is not provided, set it to nperseg
+            object.__setattr__(self, 'nfft', self.nperseg)
+
+        if not isinstance(self.window, str):
+            # If window is provided as an array, ensure it has the correct length
+            if (w := np.array(self.window).squeeze()).shape[0] != self.nperseg:
+                raise ValueError(f"Window length ({len(w)}) must match nperseg ({self.nperseg})")
+            else:
+                # If it's a valid window array, (re-)assign it as NDArray
+                object.__setattr__(self, 'window', w)
 
-        # only interpolate if needed. At least 2 points for linear interpolation are necessary
-        if idx.sum() >= 2 and (not idx.all()):
-            spec[i, idx] = np.interp(freq[idx], freq[~idx], spec[i, ~idx])
+    @property
+    def noverlap(self) -> int:
+        """
+        Calculate the number of overlapping samples based on the overlap percentage.
+        """
+        return int(self.nperseg * self.overlap)
 
-    # transpose back
-    if axis != 0:
-        spec = np.moveaxis(spec, 0, axis)
+    @property
+    def hop_length(self) -> int:
+        """
+        Calculate the hop length based on the number of FFT points and overlap.
+        """
+        return self.nperseg - self.noverlap
+
+    def get_window_correction(self, representation: ValueRepresentation) -> float:
+        win: NDArray = get_window(self.window, self.nperseg) if isinstance(self.window, str) else self.window
+        wc_db = 10 * np.log10(np.sum(win ** 2) / win.shape[0]) - 20 * np.log10(np.abs(win).mean())
+        wc_db *= -1.0 # Convert to negative value for correction
+
+        representation = representation.lower()
+        if representation == 'rms':
+            return np.power(10, wc_db / 20.0)
+        elif representation == 'pow':
+            return np.power(10, wc_db / 10.0)
+        elif representation == 'db':
+            return wc_db
+        else:
+            raise ValueError(f"Unsupported representation: {representation}")
 
-    return spec
 
-def check_spectrum_arguments(fs: int, window: str = 'flattop', nperseg: Optional[int] = None,
+def check_spectrum_arguments(fs: int, window: NDArray|str = 'flattop', nperseg: Optional[int] = None,
     noverlap: Optional[int] = None, nfft: Optional[int] = None, average: Optional[str] = None,
     fraction: Optional[int] = None,
     **kwargs) -> Dict[str, Any]:
@@ -93,8 +146,8 @@ def check_spectrum_arguments(fs: int, window: str = 'flattop', nperseg: Optional
 
     # fractional bands
     if fraction is not None:
-        Check(fraction).is_greater_than(0).or_raise(ValueError, f"Fractional bands parameter must be >= 1")
-        Check(fraction).is_an(int).or_raise(ValueError, f"Fractional bands parameter must be integer")
+        Check(fraction).is_greater_than(0).or_raise(ValueError, f"Fractional bands parameter must be >= 1, but is: {fraction}")
+        Check(fraction).is_an(int).or_raise(ValueError, f"Fractional bands parameter must be integer, but is: {fraction}")
         spec_args['fraction'] = fraction
 
     return spec_args

nonlinearity/spectrum/types.py

+"""
+Created on: 02 Jun 2025 - 14:50
+Author: jan.reimes@head-acoustics.com
+
+Type definitions and protocols for spectrum analysis functionality.
+Defines interfaces for spectrum averaging, time-varying spectra, and frequency band operations.
+Includes data classes for band definitions and protocol specifications for spectrum implementations.
+"""
+
+from typing import Protocol, Literal, Optional, Iterable, Tuple, Any, List, Self
+import numpy as np
+from scipy.signal import freqz_sos, butter
+from numpy.typing import NDArray
+from dataclasses import dataclass, field
+from enum import IntEnum
+
+CUT_OFF_ATTENUATION = -3.01  # dB
+
+class TransformFn(Protocol):
+    """Protocol for a transformation function operating on NDArray."""
+    def __call__(self, a: NDArray, *args: Any, **kwargs: Any) -> NDArray: ...
+
+AggregateFn = Literal['mean', 'median']|str|TransformFn
+Weighting = Optional[Literal['A', 'B', 'C', 'R_468', 'Z']|str]
+ValueRepresentation = Literal['rms', 'pow', 'dB']|str
+
+class Ax(IntEnum):
+    """
+    Enum for axis indices in spectrum data.
+    """
+    FREQ = 0
+    TIME = 1
+    CHANNELS = 2
+
+@dataclass(frozen=True)
+class Band:
+    """
+    Class representing a frequency band with center, lower, and upper frequencies.
+
+    Attributes:
+        center_freq (float): Center frequency of the band.
+        lower_freq (float): Lower frequency bound.
+        upper_freq (float): Upper frequency bound.
+        cutoff_attenuation (float): Attenuation at the cutoff frequencies (default: -3.01 dB).
+    """
+    center_freq: float
+    lower_freq: float
+    upper_freq: float
+    cutoff_attenuation: Optional[float] = field(default=CUT_OFF_ATTENUATION, init=True)
+
+    def __post_init__(self):
+        """
+        Validate that lower_freq < upper_freq and center_freq is within bounds.
+        Sets default cutoff_attenuation if None.
+        """
+        if self.lower_freq >= self.upper_freq:
+            raise ValueError("Lower frequency must be less than upper frequency.")
+        if self.center_freq < self.lower_freq or self.center_freq > self.upper_freq:
+            raise ValueError("Center frequency must be within the lower and upper frequencies.")
+
+        cutoff_attenuation = self.cutoff_attenuation if self.cutoff_attenuation is not None else CUT_OFF_ATTENUATION
+        object.__setattr__(self, 'cutoff_attenuation', cutoff_attenuation)
+
+    @classmethod
+    def from_center_freq_and_bandwidth(cls, center_freq: float, bandwidth: float,
+                                       cutoff_attenuation: Optional[float] = None
+                                       ) -> Self:
+        """
+        Create a Band instance from center frequency and bandwidth.
+
+        Args:
+            center_freq (float): Center frequency of the band.
+            bandwidth (float): Bandwidth of the band.
+            cutoff_attenuation (Optional[float]): Attenuation at cutoff (default: -3.01 dB).
+
+        Returns:
+            Band: New Band instance.
+        """
+
+        lower_freq = center_freq - (bandwidth / 2)
+        upper_freq = center_freq + (bandwidth / 2)
+        return cls(center_freq=center_freq, lower_freq=lower_freq, upper_freq=upper_freq,
+                   cutoff_attenuation=cutoff_attenuation
+                   )
+
+    @property
+    def bandwidth(self) -> float:
+        """
+        Calculate the bandwidth of the band.
+
+        Returns:
+            float: Bandwidth (upper_freq - lower_freq).
+        """
+
+        return self.upper_freq - self.lower_freq
+
+
+    def get_freq_resp(self, fs: int, freqs: NDArray, order: int = 24) -> NDArray:
+        """
+        Get the frequency response of the band using a bandpass filter.
+
+        Args:
+            fs (int): Sampling frequency.
+            freqs (NDArray): Frequencies at which to calculate the response.
+            order (int): Butterworth filter order (default: 24).
+
+        Returns:
+            NDArray: Magnitude response of the bandpass filter.
+        """
+        sos = butter(order, [self.lower_freq, self.upper_freq], fs=fs, btype='band', output='sos')
+        _, h = freqz_sos(sos, worN=freqs, fs=fs)
+        return np.abs(h)
+
+
+@dataclass(frozen=True)
+class FrequencyBands(List[Band]):
+    """
+    Collection of frequency bands with utilities for creation and access.
+
+    Attributes:
+        bands (list[Band]): List of Band objects.
+        title (str): Title or description of the band set.
+    """
+    bands: list[Band]
+    title: str
+
+    @classmethod
+    def from_center_freqs_and_bandwidths(cls, center_freqs: Iterable[float], bandwidths: Iterable[float],
+                                         title: str = "Frequency Bands", cutoff_attenuation: Optional[float] = None
+                                         ) -> 'FrequencyBands':
+        """
+        Create FrequencyBands from lists of center frequencies and bandwidths.
+
+        Args:
+            center_freqs (Iterable[float]): Center frequencies for each band.
+            bandwidths (Iterable[float]): Bandwidths for each band.
+            title (str): Descriptive title for the band set.
+            cutoff_attenuation (Optional[float]): Attenuation at cutoff frequencies.
+
+        Returns:
+            FrequencyBands: New FrequencyBands instance.
+
+        Raises:
+            ValueError: If center_freqs and bandwidths have different lengths.
+        """
+        if len(center_freqs) != len(bandwidths):
+            raise ValueError("Length of center_freqs and bandwidths must be the same.")
+
+        bands = [Band.from_center_freq_and_bandwidth(cf, bw, cutoff_attenuation) for cf, bw in zip(center_freqs, bandwidths)]
+        return cls(bands=bands, title=title)
+
+    @classmethod
+    def from_center_lower_upper_freqs(cls, center_freqs: Iterable[float], lower_freqs: Iterable[float],
+                                      upper_freqs: Iterable[float], title: str = "Frequency Bands",
+                                      cutoff_attenuation: Optional[float] = None
+                                      ) -> 'FrequencyBands':
+        """
+        Create FrequencyBands from center, lower, and upper frequency arrays.
+
+        Args:
+            center_freqs (Iterable[float]): Center frequencies.
+            lower_freqs (Iterable[float]): Lower frequency bounds.
+            upper_freqs (Iterable[float]): Upper frequency bounds.
+            title (str): Descriptive title for the band set.
+            cutoff_attenuation (Optional[float]): Attenuation at cutoff frequencies.
+
+        Returns:
+            FrequencyBands: New FrequencyBands instance.
+
+        Raises:
+            ValueError: If arrays have different lengths.
+        """
+        if len(center_freqs) != len(lower_freqs) or len(center_freqs) != len(upper_freqs):
+            raise ValueError("Length of center_freqs, lower_freqs, and upper_freqs must be the same.")
+
+        bands = [Band(cf, lf, uf, cutoff_attenuation) for cf, lf, uf in zip(center_freqs, lower_freqs, upper_freqs)]
+        return cls(bands=bands, title=title)
+
+    @classmethod
+    def from_band_tuples(cls, *band_tuples: Tuple[float, float, float], title: str = "Frequency Bands",
+                         cutoff_attenuation: Optional[float] = None
+                         ) -> 'FrequencyBands':
+        """
+        Create FrequencyBands from tuples of (center_freq, lower_freq, upper_freq).
+
+        Args:
+            *band_tuples: Variable number of 3-tuples (center, lower, upper frequencies).
+            title (str): Descriptive title for the band set.
+            cutoff_attenuation (Optional[float]): Attenuation at cutoff frequencies.
+
+        Returns:
+            FrequencyBands: New FrequencyBands instance.
+
+        Raises:
+            ValueError: If any tuple doesn't contain exactly 3 values.
+        """
+        if not all(len(t) == 3 for t in band_tuples):
+            raise ValueError("Each band tuple must contain three values: (center_freq, lower_freq, upper_freq).")
+
+        bands = [Band(cf, lf, uf, cutoff_attenuation) for cf, lf, uf in band_tuples]
+        return cls(bands=bands, title=title)
+
+    def __post_init__(self):
+        """Sort bands by center frequency after initialization."""
+        self.bands.sort(key=lambda band: band.center_freq)
+
+    @property
+    def center_freqs(self) -> NDArray:
+        """Get array of center frequencies for all bands."""
+        return np.array([band.center_freq for band in self.bands])
+
+    @property
+    def lower_freqs(self) -> NDArray:
+        """Get array of lower frequency bounds for all bands."""
+        return np.array([band.lower_freq for band in self.bands])
+
+    @property
+    def upper_freqs(self) -> NDArray:
+        """Get array of upper frequency bounds for all bands."""
+        return np.array([band.upper_freq for band in self.bands])
+
+    # make dataclass iterable
+    def __iter__(self) -> Iterable[Band]:
+        """Iterate over the bands in the collection."""
+        return iter(self.bands)
+
+    def __contains__(self, band: Band) -> bool:
+        """Check if a specific band is in the collection."""
+        return band in self.bands
+
+    def __getitem__(self, index: int) -> Band:
+        """Get a band by its index in the collection."""
+        return self.bands[index]
+
+    def __len__(self) -> int:
+        """Get the number of bands in the collection."""
+        return len(self.bands)
+
+# Protocol for LevelVsTime
+class ILevelVsTime(Protocol):
+    """Protocol for level-vs-time representations of audio signals."""
+
+    @property
+    def time(self) -> NDArray:
+        """Time axis values for the level data."""
+        ...
+
+    @property
+    def fs_out(self) -> float:
+        """Output sampling frequency for the level data."""
+        ...
+
+    @property
+    def axis(self) -> int:
+        """Axis along which time varies in the data arrays, default: 0 (input shape = (time, channels) or just (time,))."""
+        ...
+
+    def get_values(self, representation: ValueRepresentation = 'dB', weighting: Weighting = None) -> NDArray:
+        """
+        Get level values over time with optional weighting and representation.
+
+        Args:
+            representation (ValueRepresentation): Output format ('dB', 'rms', 'pow').
+            weighting (Weighting): Optional frequency weighting ('A', 'C', etc.).
+
+        Returns:
+            NDArray: Level values vs time.
+        """
+        ...
+
+    def get_active_frames(self, inactivity_fraction: float, weighting: Weighting = None) -> NDArray[np.bool_]:
+        """
+        Identify active frames based on inactivity threshold.
+
+        Args:
+            inactivity_fraction (float): Fraction below which frames are considered inactive.
+            weighting (Weighting): Optional frequency weighting for level calculation.
+
+        Returns:
+            NDArray[np.bool_]: Boolean mask of active frames.
+        """
+        ...
+
+# Protocols for SpectrumAvg and SpectrumVsTime.
+class ISpectrumBase(Protocol):
+    """Base protocol for all spectrum analyses (Avg and VsTime!)."""
+
+    @property
+    def _int_vr(self) -> ValueRepresentation:
+        """Internal (default) representation used for calculations."""
+        return 'pow'
+
+    @property
+    def x(self) -> NDArray:
+        """Primary input signal."""
+        ...
+
+    @property
+    def y(self) -> NDArray:
+        """Secondary input signal (for cross-spectral analysis)."""
+        ...
+
+    @property
+    def is_2d(self) -> bool:
+        """Indicates if the spectrum is always 2D (time, channels)."""
+        ...
+
+    @property
+    def fs(self) -> int:
+        """Sampling frequency of input signals."""
+        ...
+
+    @property
+    def spectrum_args(self) -> "SpectralArguments":
+        """Configuration parameters for spectral analysis."""
+        ...
+
+    def get_freqs(self, all_freqs: bool = False) -> NDArray:
+        """
+        Get frequency bins of the spectrum.
+
+        Args:
+            all_freqs (bool): If True, return all computed frequency bins.
+
+        Returns:
+            NDArray: Frequency values in Hz.
+        """
+        ...
+
+    def get_values(self, representation: ValueRepresentation = 'rms', weighting: Weighting = None, all_freqs: bool = False,
+                    value_transform: Optional[TransformFn] = None, include_window_correction: bool = True) -> NDArray:
+        """
+        Get spectrum values with optional processing and weighting.
+
+        Args:
+            representation (ValueRepresentation): Output format ('rms', 'pow', 'dB').
+            weighting (Weighting): Optional frequency weighting function.
+            all_freqs (bool): If True, return all frequency bins.
+            value_transform (Optional[TransformFn]): Optional transformation function.
+
+        Returns:
+            NDArray: Spectrum values.
+        """
+        ...
+
+class ISpectrumAvg(ISpectrumBase):
+    """Protocol for averaged spectrum functionality."""
+
+    @classmethod
+    def from_signals(cls, x: NDArray, fs: int, y: Optional[NDArray] = None, aggregate: AggregateFn = "mean", axis: int = 0, **kwargs) -> 'ISpectrumAvg':
+        """
+        Create averaged spectrum from input signals.
+
+        Args:
+            x (NDArray): Primary input signal.
+            fs (int): Sampling frequency.
+            y (Optional[NDArray]): Secondary signal for cross-spectral analysis.
+            aggregate (AggregateFn): Aggregation method ('mean', 'median', etc.).
+            axis (int): Time axis of input signal.
+            **kwargs: Additional spectral parameters.
+
+        Returns:
+            ISpectrumAvg: Averaged spectrum instance.
+        """
+        ...
+
+    def get_level(self, representation: ValueRepresentation = 'dB', weighting: Weighting = None, all_freqs: bool = False,
+                        include_window_correction: bool = True) -> float|NDArray:
+        """
+        Calculate overall level of the averaged spectrum.
+
+        Args:
+            representation (ValueRepresentation): Output format for level.
+            weighting (Weighting): Optional frequency weighting.
+            all_freqs (bool): Whether to include all frequency bins.
+
+        Returns:
+            float | NDArray: Overall level value(s).
+        """
+        ...
+
+class ISpectrumVsTime(ISpectrumBase):
+    """Protocol for time-varying spectrum (spectrogram) functionality."""
+
+    @classmethod
+    def from_signals(cls, x: NDArray, fs: int, *, y: Optional[NDArray] = None, axis: int = 0, **kwargs) -> 'ISpectrumVsTime':
+        """
+        Create time-varying spectrum from input signals.
+
+        Args:
+            x (NDArray): Primary input signal (1D or 2D).
+            fs (int): Sampling frequency.
+            y (Optional[NDArray]): Secondary signal for cross-spectral analysis.
+            axis (int): Time axis of input signal.
+            **kwargs: Additional spectral parameters.
+
+        Returns:
+            ISpectrumVsTime: Time-varying spectrum instance.
+        """
+        ...
+
+    @property
+    def time(self) -> NDArray:
+        """Time axis values for the spectrogram."""
+        ...
+
+    def get_noise_estimate(self, percentile: float = 5.0) -> ISpectrumAvg:
+        """
+        Estimate noise floor from low percentile of spectrum values.
+
+        Args:
+            percentile (float): Percentile for noise estimation (0-100).
+
+        Returns:
+            ISpectrumAvg: Estimated noise spectrum.
+        """
+        ...
+
+    def get_aggregate(self, aggregate: AggregateFn = "mean") -> ISpectrumAvg:
+        """
+        Aggregate spectrum over time using specified method.
+
+        Args:
+            aggregate (AggregateFn): Aggregation function or method name.
+
+        Returns:
+            ISpectrumAvg: Time-averaged spectrum.
+        """
+        ...
+
+    def get_level_vs_time(self) -> ILevelVsTime:
+        """
+        Calculate overall level for each time frame of the spectrum.
+
+        Returns:
+            ILevelVsTime: Level vs time representation.
+        """
+        ...
+
+    def with_active_frames(self, idx_active: NDArray[np.bool_], pad_value: bool = False) -> Self:
+        """
+        Create new instance using only specified active time frames.
+
+        Args:
+            idx_active (NDArray[np.bool_]): Boolean mask indicating active frames.
+            pad_value (bool): Whether to pad inactive frames with specified values.
+
+        Returns:
+            ISpectrumVsTime: New instance with active frames.
+        """
+        ...
+
+
+if __name__ == "__main__":
+    pass
\ No newline at end of file

nonlinearity/spectrum/weighting.py

+"""
+Created on: 04 Jun 2025 - 10:45
+Author: jan.reimes@head-acoustics.com
+
+This module provides frequency weighting functions for acoustic measurements.
+Supports standard weightings (A, C, Z) as well as deprecated (B, D) and
+specialized weightings (ITU-R 468) commonly used in audio and acoustics.
+
+"""
+import numpy as np
+from numpy.typing import NDArray
+from typing import Literal
+import acoustics.standards.iec_61672_1_2013 as iec_61672_1_2013
+
+from .types import Weighting, ValueRepresentation
+from .misc import get_value_representation
+
+
+def _deprecated_weighting_function(freq: NDArray|float, w: Literal['B', 'D']|str) -> NDArray|float:
+    """
+    Calculate B-/D-weighting in dB for a given frequency f (Hz).
+    References:
+        https://en.wikipedia.org/wiki/A-weighting#Function_realisation_of_some_common_weightings
+        https://www.sps.tue.nl/rmaarts/RMA_papers/aar92a.pdf
+    """
+
+    f2 = np.square(freq)
+
+    if w.lower() == 'b':
+        cb = 9.8767069950664E-1
+        cc = 6.6709544848173E-9
+
+        ps1b = np.square(p1b := 158.5)
+        ps1c = np.square(p1c := 20.6)
+        ps2c = np.square(p2c := 12200.0)
+
+        cw = f2 / (f2 + ps1c) / (f2 + ps2c) / cc
+        return 20*np.log10(cw * freq  / np.sqrt((f2 + ps1b) / cb))
+
+    elif w.lower() == 'd':
+        cd = 6.8966888496476E-5
+
+        ps1d = np.square(p1d := 282.7)
+        ps2d = np.square(p2d := 1160.0)
+        ps3d1 = np.square(p3d1 := 1712.0)
+        ps3d2 = np.square(p3d2 := 2628.0)
+        zs1d1 = np.square(z1d1 := 519.8)
+        zs1d2 = np.square(z1d2 := 876.2)
+
+        h = 1/(np.square(ps3d1 + ps3d2 - f2) + (4*f2*ps3d1))
+        h *= (np.square(zs1d1 + zs1d2 - f2) + (4*f2*zs1d1))
+        return 20*np.log10(freq  * np.sqrt(h / (f2 + ps1d) / (f2 + ps2d)) / cd)
+
+    else:
+        raise ValueError(f"Unknown weighting function: {w}")
+
+
+def get_weighting(freq: NDArray, weighting: Weighting, representation: ValueRepresentation = 'rms') -> NDArray:
+    """
+    Calculates the frequency weighting curve for a given set of frequencies and weighting standard.
+
+    Parameters
+    ----------
+    freq : NDArray
+        Array of frequency values (in Hz) for which the weighting should be calculated.
+    weighting : Weighting
+        The type of frequency weighting to apply. Supported values are:
+        - 'a', 'c': Standard A/C-weighting (via acoustics module).
+        - 'b', 'd': Deprecated B/D-weighting.
+        - 'itu_r_468': ITU-R 468 weighting.
+        - 'z' or None: No weighting (returns zeros).
+    representation : ValueRepresentation, optional
+        The representation of the output values. Default is 'rms'.
+
+    Returns
+    -------
+    NDArray
+        The weighting values for the given frequencies, in the specified representation (default: dB).
+
+    Raises
+    ------
+    ValueError
+        If an unknown weighting function is specified.
+    """
+
+    if weighting is None or weighting.lower() == 'z':
+        # no weighting
+        w = np.zeros_like(freq)
+    elif weighting.lower() in ['a', 'c']:
+        # A/C-weighting via acoustics module
+        fn_w = getattr(iec_61672_1_2013, f'weighting_function_{weighting.lower()}')
+        w = fn_w(freq)
+
+    elif weighting.lower() in ['b', 'd']:
+        # deprecated B/D-weighting
+        w = _deprecated_weighting_function(freq, weighting)
+
+    elif weighting.lower() == 'itu_r_468':
+        # ITU-R 468 weighting
+        f = {i+1: np.power(freq, i+1) for i in range(6)}
+        h1 = (-4.7373389813783836e-24 * f[6]) + (2.0438283336061252e-15 * f[4]) - (1.363894795463638e-07 * f[2]) + 1
+        h2 = (1.3066122574128241e-19 * f[5]) - (2.1181508875186556e-11 * f[3]) + (0.0005559488023498643 * f[1])
+        r_itu = (0.0001246332637532143 * f[1]) / np.sqrt(h1 * h1 + h2 * h2)
+
+        w = 18.246265068039158 + 20.0 * np.log10(r_itu)
+
+    else:
+        raise ValueError(f"Unknown weighting function: {weighting}")
+
+    return get_value_representation(w, representation, 'dB')
+
+
+if __name__ == "__main__":
+    pass
\ No newline at end of file

nonlinearity/utils/__init__.py

@@ -5,5 +5,62 @@ Created on May 29 2024 16:34
 @author: Jan.Reimes
 """
 
+from pathlib import Path
+from typing import AnyStr, Iterable, List, Tuple
+
+PathLike = str | Path
+FileInput = AnyStr | Iterable[AnyStr]
+
+
+def get_file_pairs(measured: FileInput, ref: FileInput) -> Tuple[Tuple[str, str]]:
+    # allow:
+    # - measured / ref are strings (1->1)
+    # - one of measured / ref is a list of strings, the other one is a string (1->N)
+    # - both measured / ref are lists of strings (N->N)
+
+    if isinstance(measured, PathLike):
+        measured: List[str] = [str(measured)]
+    elif not isinstance(measured, Iterable):
+        raise ValueError("Input 'measured' must be a string or an iterable of strings.")
+
+    lm = len(measured)
+
+    if isinstance(ref, PathLike):
+        ref: List[str] = [str(ref)]
+    elif not isinstance(ref, Iterable):
+        raise ValueError("Input 'ref' must be a string or an iterable of strings.")
+
+    lr = len(ref)
+
+    if lm == lr:
+        pairs = [(measured[i], ref[i]) for i in range(lm)]
+    elif (min(lm, lr) == 1) and (max(lm, lr) > 1):
+        # one of them has only one item, so we can pair it with all items of the other
+        if lm == 1:
+            pairs = [(measured[0], ref[i]) for i in range(lr)]
+        else:
+            pairs = [(measured[i], ref[0]) for i in range(lm)]
+    else:
+        raise ValueError(
+            "Inputs 'measured' and 'reference' files must have the same number of entries (or one of them has only one item)."
+        )
+
+    return (*pairs,)
+
+
 if __name__ == "__main__":
-    pass
+    pairs1 = get_file_pairs("measured.wav", "ref.wav")
+    pairs2 = get_file_pairs(["measured.wav"], ["ref.wav"])
+    print(pairs1 == pairs2)  # should be True
+
+    pairs1 = get_file_pairs(["measured1.wav", "measured2.wav"], "ref.wav")
+    pairs2 = get_file_pairs(["measured1.wav", "measured2.wav"], ["ref.wav"])
+    print(pairs1 == pairs2)  # should be True
+
+    pairs1 = get_file_pairs("measured.wav", ["ref1.wav", "ref2.wav"])
+    pairs2 = get_file_pairs(["measured.wav"], ["ref1.wav", "ref2.wav"])
+    print(pairs1 == pairs2)  # should be True
+
+    pairs1 = get_file_pairs("measured.wav", "ref.wav")
+    pairs2 = get_file_pairs("measured.wav", ["ref.wav"])
+    print(pairs1 == pairs2)  # should be True

nonlinearity/utils/p56/asl.py

@@ -4,22 +4,28 @@ Created on Jun 23 2021 13:04
 
 @author: Jan.Reimes
 
-Warning: This implementation might not 100% comply with ITU-T P.56!
+Warning: This implementation might not 100% comply with reference implementation in ITU-T P.56!
+
+Implementation adapted from:
+https://github.com/nii-yamagishilab/NELE-GAN/blob/master/asl_P56.py
+
 """
 
-from typing import Tuple
+from typing import Tuple, Optional
+from soxr import resample
 import numpy as np
 from numpy.typing import NDArray
-# from numba import jit
-from scipy.signal import lfilter
+from scipy.signal import lfilter, sosfilt, butter
 
 from .prefilter import P56Prefilter, PrefilterP56, get_filter, apply_filters
 
 class ASLException(Exception):
     pass
 
+NDArray32f = NDArray[np.float32]
+NDArray64i = NDArray[np.int64]
+
 
-#@jit(nopython=True)
 def __bin_interp(upcount: int, lwcount: int, upthr: float, lwthr: float, Margin: float, tol: float) -> Tuple[float, float]:
     tol = np.abs(tol)
 
@@ -64,7 +70,36 @@ def __bin_interp(upcount: int, lwcount: int, upthr: float, lwthr: float, Margin:
 
         return asl_ms_log, cc
 
-def calculate_p56_asl(x: NDArray, fs: int, nbits: int = 16, M: float = 15.9, H: float = 0.2, T: float = 0.03) -> Tuple[float, float]:
+def __get_activity(q: NDArray32f, c: NDArray32f, I: int) -> Tuple[NDArray64i, NDArray64i]:
+    # activity counter for each level threshold
+    #thres_no = np.int64(c.shape[0])
+    thres_no = int(c.shape[0])
+    I = np.int64(I)
+
+    a: NDArray64i = np.zeros(thres_no, dtype=np.int64)
+
+    # hangover counter for each level threshold
+    hang: NDArray64i = np.zeros(thres_no, dtype=np.int64) + I
+
+    for k in np.arange(int(x_len := q.shape[0]), dtype=int):
+        for j in np.arange(thres_no):
+            if q[k] >= c[j]:
+                # q is active -> reset hangover counter
+                a[j] = a[j] + 1
+                hang[j] = 0
+            else:
+                if hang[j] < I:
+                    # q still counted as active (in hangover time) -> increment update hangover counter
+                    a[j] = a[j] + 1
+
+                hang[j] = hang[j] + 1
+                #else:
+                #    break
+
+    return a, hang
+
+def calculate_p56_asl(x: NDArray, fs: int, nbits: int = 16, M: float = 15.9, H: float = 0.2, T: float = 0.03,
+                      **kwargs) -> Tuple[float, float]:
     '''
     This implements ITU P.56 method B.
     Usage:  asl_P56(x, fs, nbits)
@@ -73,49 +108,48 @@ def calculate_p56_asl(x: NDArray, fs: int, nbits: int = 16, M: float = 15.9, H:
         nbits         - the number of bits (default: 16)
         M             - margin in dB of the difference between threshold and active speech level (default: 15.9)
         H             - hangover time in seconds (default: 0.2)
-        T             - time constant of smoothing, in seconds (default:0.3
+        T             - time constant of smoothing, in seconds (default:0.3)
+        f0_hp         - optional: high-pass filter to be applied on x, to e.g., remove low-frequency noise
+        order_hp      - order of IIR high-pass filter
+
     Example call:
         asl_rms, asl, c0 = asl_P56(x, fs, nbits)
+
     References:
     [1] ITU-T (1993). Objective measurement of active speech level. ITU-T
         Recommendation P. 56
     Python implementation from MATLAB: Rui Cheng
     '''
 
-    thres_no = nbits - 1  # number of thresholds, for 16 bit, it's 15
-    eps = 2.2204e-16
+    # check/prepare input signal x
+    x = x.astype('float32')
+    x -= x.mean()
+    f0_hp: float = kwargs.get('f0_hp', 20.0)
+    order_hp: int = kwargs.get('order_hp', 16)
+    if f0_hp is not None:
+        sos = butter(order_hp, f0_hp, 'high', output='sos', fs=fs)
+        x = sosfilt(sos, x)
 
+    eps = np.finfo(np.float32).eps # 2.2204e-16
+
+    # analysis parameter
+    thres_no = nbits - 1  # number of thresholds, for 16 bit, it's 15
     I = int(np.ceil(fs * H))  # hangover in samples
-    g = np.exp(-1 / (fs * T))  # smoothing factor in enevlop detection
-    c = np.array([pow(2, i) for i in range(-thres_no, thres_no - nbits + 1)])
-    # vector with thresholds from one quantizing level up to half the maximum code, at a step of 2, in the case of 16bit samples, from 2^-15 to 0.5
-    a = np.zeros(thres_no, dtype=int) # activity counter for each level threshold
-    hang = I + np.zeros(thres_no, dtype=int)  # % hangover counter for each level threshold
+    g = np.exp(-1 / (fs * T))  # smoothing factor in envelop detection
 
-    sq = np.sum(np.power(x, 2))  # long-term level square energy of x
+    # vector with thresholds from one quantization level up to half the maximum code, at a step of 2, in the case of 16bit samples, from 2^-15 to 0.5
+    c = np.array([pow(2, i) for i in range(-thres_no, thres_no - nbits + 1)], dtype=np.float32)
+
+    # long-term level square energy of x
+    sq = np.sum(np.power(x, 2))
     x_len = len(x)  # length of x
 
     # use a 2nd order IIR filter to detect the envelope q
     x_abs = np.abs(x)
-    p = lfilter([1 - g], [1, -g], x_abs)
-    q = lfilter([1 - g], [1, -g], p)
-
-    #@jit(nopython=True)
-    def __get_activity(q, c, a, hang, thres_no, x_len):
-        for k in range(x_len):
-            for j in range(thres_no):
-                if q[k] >= c[j]:
-                    a[j] = a[j] + 1
-                    hang[j] = 0
-                elif hang[j] < I:
-                    a[j] = a[j] + 1
-                    hang[j] = hang[j] + 1
-                else:
-                    break
-
-        return a, hang
+    p = lfilter([1 - g], [1, -g], x_abs).astype('float32')
+    q = lfilter([1 - g], [1, -g], p).astype('float32')
 
-    a, hang = __get_activity(q, c, a, hang, thres_no, x_len)
+    a, hang = __get_activity(q, c, I)
 
     # default result values
     activity = 0
@@ -157,6 +191,12 @@ def calculate_p56_asl(x: NDArray, fs: int, nbits: int = 16, M: float = 15.9, H:
     return asl_dB, activity
 
 def calculate_p56_asl_ex(x: NDArray, fs: int, pre_filter: PrefilterP56 = 'NoFilter', min_amplitude: float = 0.1, max_amplitude: float = 1.0, **kwargs) -> Tuple[float, float]:
+    # check downsampling
+    fs_down = kwargs.get('fs_down', 16000)
+    if fs_down < fs:
+        x = resample(x, fs, fs_down)
+        fs = fs_down
+
     # call calculateP56ASL() with additional pre-filter and range check of signal:
     x = np.array(x)
     max_abs_value = np.abs(x).max()
@@ -166,6 +206,8 @@ def calculate_p56_asl_ex(x: NDArray, fs: int, pre_filter: PrefilterP56 = 'NoFilt
     if pre_filter != P56Prefilter.NoFilter:
         coeffs = get_filter(pre_filter, fs)
         y = apply_filters(x, coeffs)
+        # pre-filter already applied -> disable hp20 in main calculation explicitly
+        kwargs['f0_hp'] = None
     else:
         y = x
 

nonlinearity/utils/resolver.py

+from pathlib import Path
+from functools import partial
+from typing import List
+from omegaconf import OmegaConf
+
+
+def _pathlib_glob(data_dir: str, *glob_patterns: str, recursive: bool = True) -> List[str]:
+    """Use Pathlib (r)glob to get a list of filenames"""
+
+    data_dir_path = Path(data_dir)
+    glob_cmd = data_dir_path.rglob if recursive else data_dir_path.glob
+    file_paths: List[Path] = []
+    for glob_pattern in glob_patterns:
+        file_paths += list(glob_cmd(glob_pattern))
+
+    filenames: List[str] = [str(p) for p in file_paths]
+    return filenames
+
+path_glob = partial(_pathlib_glob, recursive=False)
+path_rglob = partial(_pathlib_glob, recursive=True)
+
+
+def register_all_resolvers():
+    OmegaConf.register_new_resolver("path_glob", path_glob)
+    OmegaConf.register_new_resolver("path_rglob", path_rglob)
+
+
+if __name__ == "__main__":
+    # usage examples:
+    register_all_resolvers()
+
+    yaml_data = r"""
+    paths:
+        data_dir: d:\\SVN\\Projects\\ts104063-nonlinearity-measure\\data\\TR-41\\All Recordings
+
+    file_names1: ${pathlib_rglob:${paths.data_dir}, '*.wav'}
+    file_names2: ${pathlib_rglob:${paths.data_dir}, 'DUT A/**/*.wav'}
+
+    """
+    cfg = OmegaConf.create(yaml_data)
+    print(len(cfg.file_names1))
+    print(len(cfg.file_names2))
+
+    pass

pyproject.toml

@@ -2,39 +2,69 @@
 name = "Nonlinearity-Measure"
 version = "0.1.0"
 description = "Code for non-linearity measure for upcoming ETSI TS 104 063"
-requires-python = ">=3.11,<3.13"
+requires-python = ">=3.12,<=3.13"
 dependencies = [
     "colorednoise>=2.2.0",
     "ensure>=1.0.4",
     "hydra-core>=1.3.2",
     "hydra_colorlog",
-    "matplotlib>=3.9.3",
     "numexpr>=2.10.2",
     "numpy>=2.1.3",
     "omegaconf>=2.3.0",
-    "pandas[excel, hdf5]>=2.2.3",
+    "pandas>=2.2.3",
     "rootutils>=1.0.7",
     "scipy>=1.14.1",
     "soundfile>=0.12.1",
     "soxr>=0.5.0.post1",
-    "wxPython",
-    "openpyxl",
     "hydra-list-sweeper>=1.0.1",
     "hydra-callbacks>=0.6.1",
-    "h5py>=3.12.1",
     "hydra-joblib-launcher>=1.2.0",
+    "python-calamine>=0.3.1",
+    "xlsxwriter>=3.2.2",
+    "dependency-groups>=1.3.0",
+    "numba>=0.61.2",
+    "acoustics>=0.2.6",
+    "tqdm>=4.67.1",
+    "tqdm-joblib>=0.0.4",
+    "psutil>=7.0.0",
+]
+
+[project.optional-dependencies]
+debug = [
     "locket>=1.0.0",
+    "h5py>=3.13.0",
+    "matplotlib>=3.10.1",
+    "seaborn>=0.13.2",
+    "mplcursors>=0.6",
 ]
 
 [dependency-groups]
 dev = [
-    "acoustics>=0.2.6",
     "ddt>=1.7.2",
+    "hydra-yaml-lsp>=0.1.1",
+    "nuitka>=2.7",
     "psutil>=6.1.1",
     "pytest>=8.3.4",
+    "pytest-cov>=6.2.1",
     "pyyaml>=6.0.2",
-    "seaborn>=0.13.2",
 ]
 
+[tool.pytest]
+log_cli=true
+log_level="INFO"
+
 [tool.pytest.ini_options]
 pythonpath = "nonlinearity_measure"
+addopts = [
+    "--import-mode=importlib",
+]
+
+[tool.uv]
+link-mode = "symlink"
+
+[tool.pyright]
+typeCheckingMode = "off"
+ignore = ["*"]
+
+[tool.ruff]
+line-length = 188

pytest.ini

-[pytest]
-log_cli=true
-log_level=INFO
\ No newline at end of file

tests/__init__.py

@@ -6,6 +6,7 @@ Created on May 29 2024 15:58
 """
 
 from pathlib import Path
+from typing import Optional, AnyStr
 from rootutils import setup_root
 
 test_dir = Path(__file__).parent
@@ -13,11 +14,17 @@ _output_dir = test_dir / 'output'
 
 setup_root(__file__, pythonpath=True)
 
-def get_output_dir(create_if_not_existent: bool = True) -> Path:
+PathLike = AnyStr | Path
+
+def get_output_dir(sub_dir: Optional[PathLike] = None, create_if_not_existent: bool = True) -> Path:
+    the_output_dir = _output_dir
+    if sub_dir is not None:
+        the_output_dir = _output_dir / sub_dir
+
     if create_if_not_existent:
-        _output_dir.mkdir(exist_ok=True, parents=True)
+        the_output_dir.mkdir(exist_ok=True, parents=True)
 
-    return _output_dir
+    return the_output_dir
 
 if __name__ == "__main__":
     pass

tests/test_active_frames.py

@@ -8,8 +8,9 @@ from itertools import groupby
 import pandas as pd
 from scipy.signal import get_window
 
+from nonlinearity import to_db
 from nonlinearity.utils.p56.asl import calculate_p56_asl_ex
-from nonlinearity.spectrum import csd_vs_time
+from nonlinearity.spectrum import FftSpectrumVsTime
 
 from tests import get_output_dir
 from tests.testfiles.file_dl import get_p501_file
@@ -54,7 +55,9 @@ class ActiveFramesTestCase(unittest.TestCase):
         inactivity = 1 - act
 
         # use auto power spectra vs time to calculate level vs time
-        f, t, pxx = csd_vs_time(x, x, fs=fs, **spec_args)
+        spec = FftSpectrumVsTime.from_signals(x, y=x, fs=fs, **spec_args)
+        pxx = spec.get_values('pow')
+        f, t = spec.get_freqs(), spec.time
 
         idx_f = np.ones_like(f, dtype=bool)
         if fmin is not None:
@@ -63,7 +66,7 @@ class ActiveFramesTestCase(unittest.TestCase):
         if fmax is not None:
             idx_f &= f <= fmax
 
-        lxx_db = 10*np.log10(np.maximum(np.sum(np.abs(pxx[idx_f, :]), axis=0), 1e-10)) + wcf
+        lxx_db = to_db(np.sum(np.abs(pxx[idx_f, :]), axis=0), is_rms=False) + wcf
 
         # get inactivity threshold
         activity_threshold_db = np.percentile(lxx_db, 100 * inactivity)

tests/test_base.py

@@ -5,7 +5,7 @@ Created on Okt 03 2024 02:32
 @author: Jan.Reimes
 """
 import unittest
-from omegaconf import OmegaConf, DictConfig
+from omegaconf import DictConfig
 from pathlib import Path
 from hydra.core.hydra_config import HydraConfig
 from hydra import initialize, compose
@@ -14,14 +14,14 @@ from typing import Union, Iterable, Tuple, Optional, List, Any, Dict, Generator
 import numpy as np
 from numpy.typing import NDArray
 from soxr import resample
-from contextlib import contextmanager, AbstractContextManager
+from contextlib import contextmanager
 
 from nonlinearity.utils.filter import lfilter_block
+from nonlinearity.spectrum import FftSpectrumAvg
 from tools.harmonics import create_harmonics
 
 from tests.testfiles.file_dl import get_g191_file
 
-MIN_DB = -120.0
 
 @contextmanager
 def get_config(overrides: Optional[List[str]] = None, no_custom_config: bool = True, config_name: str = "config.yaml",
@@ -44,30 +44,31 @@ class BaseTestCase(unittest.TestCase):
     @classmethod
     def setUpClass(cls):
         # auto setup for making all submodules relative to this available file
-        cls.root_dir = setup_root(__file__, cwd=True)
+        cls.root_dir = setup_root(__file__, cwd=True, pythonpath=True)
         cls.test_dir = Path(__file__).parent
 
     @classmethod
-    def _distort_signal(cls, source: NDArray, fs: int = None, delay=0.04, snr: float = 50.0,
+    def _distort_signal(cls, source: NDArray, fs: int, delay=0.04, snr: float = 50.0,
                         nbr_channels: int = 2,  #up_down_sampling_factor: int = 3,
                         harmonics: Union[int, Iterable[int]] = 4,
                         harmonic_level_offsets_db: float = -10.0,
-                        overall_level_harmonics_db: float = None,
+                        overall_level_harmonics_db: Optional[float] = None,
                         max_harmonic_freq: Optional[float] = None,
                         min_harmonic_freq: Optional[float] = None,
                         rir_filter_id: str = "LEABP01") -> Tuple[NDArray, NDArray]:
 
-        fs = cls.cfg.fs if fs is None else fs
-        level_db = 20 * np.log10(np.abs(source).max() / np.sqrt(2))
-
-        # upsample by factor up_down_sampling_factor (avoid alias components of distortion )
-        #if up_down_sampling_factor > 1:
-        #    s = resample(source, fs, up_down_sampling_factor*fs, 'VHQ')
+        if (nbr_channels > 2) and (rir_filter_id is not None):
+            raise ValueError("Only 2 channels are supported for usage with rir_filter.")
 
         d = round(delay * fs)
-        harmonic_components = np.zeros((source.shape[0] + d, nbr_channels))
-        processed = harmonic_components.copy()
+
+        processed = np.zeros((source.shape[0] + d, nbr_channels))
+        harmonic_components = np.zeros_like(processed)
+
         for i in range(nbr_channels):
+            processed[d:, i] = source.copy()
+
+            # apply harmonic distortions
             harmonic_components[d:, i] = create_harmonics(source, fs, harmonics,
                                                  max_harmonic_freq=max_harmonic_freq,
                                                  min_harmonic_freq=min_harmonic_freq,
@@ -76,18 +77,10 @@ class BaseTestCase(unittest.TestCase):
                                                  harmonics_only=True,
                                                  )
 
-            processed[d:, i] = harmonic_components[d:, i] + source
-
-        # downsample by factor up_down_sampling_factor (avoid alias components of distortion )
-        #if up_down_sampling_factor > 1:
-        #    sos = sig.butter(8, 0.95*(fs/2), 'low', output='sos', fs=up_down_sampling_factor*fs)
-        #    processed = sig.sosfiltfilt(sos, processed, axis=0)
-        #    processed = resample(processed, up_down_sampling_factor*fs, fs, 'VHQ')
-
             if rir_filter_id is not None:
                 # apply RIR from G.191
-            for i, ch in enumerate(['L', 'R']):
-                if i < nbr_channels:
+                ch = 'LR'[i]
+
                 rir_file = get_g191_file(fr'little_endian\{rir_filter_id}.{ch}.IR32')
                 with rir_file.open('rb') as f:
                     rir = np.fromfile(f, dtype=np.float32)
@@ -95,15 +88,18 @@ class BaseTestCase(unittest.TestCase):
                 if fs != 32000:
                     rir = resample(rir, 32000, fs)
 
-                    processed[:, i] = lfilter_block(rir, np.array([1.0]), processed[:, i], compensate_fir_delay=False)
-                    harmonic_components[:, i] = lfilter_block(rir, np.array([1.0]), harmonic_components[:, i], compensate_fir_delay=False)
+                for x in [processed, harmonic_components]:
+                    x[:, i] = lfilter_block(rir, np.array([1.0]), x[:, i], compensate_fir_delay=False)
+
+            #processed += harmonic_components
 
         # idle noise
+        level_db = FftSpectrumAvg.from_signals(processed, fs).get_level('dB')
+        #level_db = calculate_level(processed, fs)
         noise_level_db = level_db - snr
         noise = np.random.randn(*processed.shape) * np.power(10, noise_level_db / 20)
-        processed += noise
 
-        return processed, harmonic_components
+        return processed + harmonic_components + noise, harmonic_components
 
     @staticmethod
     def _pad_same_length(*signals: NDArray) -> Tuple[NDArray, ...]:
@@ -121,102 +117,5 @@ class BaseTestCase(unittest.TestCase):
         return tuple(padded)
 
 
-    @staticmethod
-    def fft_to_octave_bands(fft_data, fs, fraction=1, is_transfer_function=False, base_freq=1000.0,
-                            clip_output_db=-100.0):
-        """
-        Transforms FFT frequency response into fractional octave bands with linear interpolation,
-        handling transfer functions and optional output clipping in dB.
-
-        Parameters:
-        - fft_data: Array of FFT magnitudes or power spectrum (assumed to be real FFT, i.e., output of np.fft.rfft)
-        - fs: Sampling frequency in Hz
-        - fraction: Fractional octave (e.g., 1 for full octave, 3 for 1/3 octave)
-        - is_transfer_function: Boolean flag. If True, applies averaging for transfer function frequency responses.
-        - base_freq: The reference center frequency for octave bands (default is 1000 Hz)
-        - clip_output_db: Minimum dB level for output clipping (default is -100 dB)
-
-        Returns:
-        - octave_band_levels: Array with the magnitude/power in each octave band (clipped at clip_output_db)
-        - band_center_freqs: Corresponding center frequencies for the octave bands
-        """
-
-        # Ensure we work with magnitudes (handle complex fft_data)
-        fft_data = np.abs(fft_data)
-
-        # Derive n_fft from fft_data size
-        n_fft = (len(fft_data) - 1) * 2
-
-        # Frequency vector for FFT bins
-        freqs = np.fft.rfftfreq(n_fft, 1 / fs)
-
-        # Calculate clipping level in linear domain from dB
-        clip_level_linear = 10 ** (clip_output_db / 20.0)  # dB to linear conversion (assuming magnitudes)
-
-        # Define constants for fractional octave bands
-        n_bands = np.ceil(np.log2(freqs[-1] / base_freq) * fraction)
-
-        # Calculate center frequencies of the octave bands
-        band_center_freqs = base_freq * 2 ** (np.arange(-n_bands, n_bands + 1) / fraction)
-
-        # Initialize an array to store the magnitude/power in each octave band
-        octave_band_levels = np.zeros(len(band_center_freqs))
-
-        # Loop over each octave band
-        f_center: float
-        for i, f_center in enumerate(band_center_freqs):
-            # Band edges
-            f_low = f_center / (2 ** (1 / (2 * fraction)))
-            f_high = f_center * (2 ** (1 / (2 * fraction)))
-
-            # Find the FFT bins that fall within this frequency band
-            idx_band = np.where((freqs >= f_low) & (freqs <= f_high))[0]
-
-            if len(idx_band) > 0:
-                if is_transfer_function:
-                    # Average the magnitudes for transfer functions
-                    octave_band_levels[i] = np.mean(fft_data[idx_band])
-                else:
-                    # Sum the magnitudes/power (use squared magnitudes for power)
-                    octave_band_levels[i] = np.sqrt(np.sum(fft_data[idx_band] ** 2))
-            else:
-                # Linear interpolation when no FFT bin falls within the band
-                if f_low < freqs[0]:
-                    # If the lower frequency bound is below the first FFT bin, set to zero
-                    octave_band_levels[i] = 0
-                elif f_high > freqs[-1]:
-                    # If the upper frequency bound is beyond the highest FFT bin, set to zero
-                    octave_band_levels[i] = 0
-                else:
-                    # Perform linear interpolation
-                    lower_bin = np.searchsorted(freqs, f_low, side='right') - 1
-                    upper_bin = np.searchsorted(freqs, f_high, side='left')
-
-                    if lower_bin >= 0 and upper_bin < len(freqs):
-                        # Interpolate between the nearest bins
-                        interpolated_value = np.interp(
-                            [f_center],
-                            [freqs[lower_bin], freqs[upper_bin]],
-                            [fft_data[lower_bin], fft_data[upper_bin]]
-                        )[0]
-
-                        if is_transfer_function:
-                            # Use interpolated value for transfer functions
-                            octave_band_levels[i] = interpolated_value
-                        else:
-                            # Sum the interpolated value for non-transfer functions
-                            octave_band_levels[i] = np.sqrt(interpolated_value ** 2)
-
-        # Apply clipping in the linear domain (convert clip_output_db to linear)
-        octave_band_levels = np.maximum(octave_band_levels, clip_level_linear)
-
-        return octave_band_levels, band_center_freqs
-
-    @staticmethod
-    def fr_to_db(fr, min_db=MIN_DB, is_rms: bool = True):
-        fr = np.abs(fr) # .squeeze()
-        fr = np.maximum(fr, np.power(10, min_db / (db_factor := 20.0 if is_rms else 10.0)))
-        return db_factor * np.log10(fr)
-
 if __name__ == "__main__":
     pass

tests/test_calculate_nlm_main.py

@@ -6,15 +6,18 @@ Created on Sept 05 2024 18:22
 """
 import unittest
 from pathlib import Path
-import numpy as np
 from ddt import ddt, data
 
-from nonlinearity import NLMResult, MIN_LEVEL_LIN
+from nonlinearity import NLMResult
 from nonlinearity.__main__ import main
 
+from tools.plot_spectra import NLMPlot
+
 from tests.test_base import BaseTestCase, get_config
+from tests import get_output_dir
 from tests.testfiles.file_dl import get_p501_file
 
+
 @ddt
 class NlmCalculateTestCase(BaseTestCase):
 
@@ -26,21 +29,16 @@ class NlmCalculateTestCase(BaseTestCase):
         overrides = [f'file_measured="{deg}"', f'file_ref="{ref}"', f"use_activity_threshold={use_activity_threshold}"]
 
         # run with config
-        with get_config(overrides=overrides) as cfg:
+        with get_config(overrides=overrides, no_custom_config=False) as cfg:
             res: NLMResult = main(cfg)
 
             # plot results
-        import matplotlib.pyplot as plt
-        levels = res.levels.asdict()
-        for name, spectrum in res.spectra.asdict().items():
-            spectrum = np.maximum(spectrum[:,0], MIN_LEVEL_LIN)
-            
-            plt.semilogx(res.freq, 20 * np.log10(spectrum.squeeze()), label=f"{name} [{levels[name][0]:.1f} dB]")
+            nlm_plot = NLMPlot(res, title=f"NLM Test: {deg.name}")
+            nlm_plot.plot(plot_nlm=True, plot_fr=True, plot_components=True, plot_cxy=True,
+                        show_blocking=True,
+                        save_path=get_output_dir(sub_dir="plots")
+                        )
 
-        plt.grid(True)
-        plt.title(f"NLM = {res.nlm.squeeze():.1f} dB (use_activity_threshold={use_activity_threshold})")
-        plt.legend()
-        plt.show()
 
 
 if __name__ == '__main__':

tests/test_ir_estimate.py

+import unittest
+from pathlib import Path
+from numpy.typing import NDArray
+import numpy as np
+import soundfile as sf
+import matplotlib.pyplot as plt
+from scipy.signal import butter, sosfilt, freqz_sos
+from ddt import ddt, data
+
+from nonlinearity.estimate_ir import estimate_ir_csd
+from nonlinearity.utils.p56 import calculate_p56_asl_ex
+from nonlinearity.spectrum.types import FrequencyBands
+from nonlinearity.spectrum import get_value_representation
+from nonlinearity.spectrum.bands import get_bands_iec_61260_1_2014
+
+from tests import get_output_dir
+from tests.testfiles.file_dl import get_p501_file
+
+@ddt
+class IrEstimateTestCase(unittest.TestCase):
+
+    @classmethod
+    def setUpClass(cls):
+        # Set the output directory for plots
+        cls.output_dir = get_output_dir("ir_estimate")
+
+        # Load the reference audio file
+        #ref: Path = get_p501_file('English_FB_clause_7.3/FB_male_female_single-talk_seq.wav')
+        ref: Path = get_p501_file('P501_D_AM_fm_FB_48k.wav')
+        cls.s, cls.fs = sf.read(ref)
+        cls.asl, cls.act = calculate_p56_asl_ex(cls.s, cls.fs, 'SWB')
+
+    def _get_noise_signal(self, snr: float = 30.0) -> NDArray:
+        return 10**((self.asl - snr)/20) * np.random.randn(*self.s.shape)
+
+    @data(True, False)
+    def test_ir_estimate(self, use_activity: bool, snr: float = 15.0, nfft: int = 16384, order: int = 4, ntaps: int = 2*4096) -> None:
+        inactivity_fraction = (1-self.act) if use_activity else None
+
+        # target fft frequency vector
+        freqs = np.fft.rfftfreq(nfft, d=1/self.fs)
+        s = self.s + self._get_noise_signal(snr)
+        fs = self.fs
+
+        # get frequency bands
+        bands: FrequencyBands = get_bands_iec_61260_1_2014(fraction=3, fmin=100.0, fmax=10_000.0)
+        fig, ax = plt.subplots(1,1, figsize=(14, 6))
+
+        for i, band in enumerate(bands):
+            with self.subTest(band=f"[{(i+1):02d}] {band.lower_freq}-{band.upper_freq} Hz"):
+                # get filter coeffs and apply to source signal
+                sos = butter(order, [band.lower_freq, band.upper_freq], fs=fs, btype='band', output='sos')
+                y = sosfilt(sos, s)
+
+                # real frequency response of the filter
+                _, h = freqz_sos(sos, worN=freqs, fs=fs)
+
+                # Estimate the impulse response using cross-spectral density
+                ir_res = estimate_ir_csd(y, s, fs, nperseg=nfft, noverlap=round(0.75*nfft), ntaps=ntaps, inactivity_fraction=inactivity_fraction)
+
+                # Check the results
+                self.assertIsNotNone(ir_res, "Impulse response estimation failed.")
+
+                # delta in the frequency range of each band must be below threshold
+                h_db = get_value_representation(h, 'dB', 'rms')
+                est_db = get_value_representation(ir_res.freq_resp, 'dB', 'rms')
+                idx_f1 = np.argmin(np.abs(freqs - band.lower_freq))
+                idx_f2 = np.argmin(np.abs(freqs - band.upper_freq))+1
+                delta = np.around(np.abs(h_db - est_db)[idx_f1:idx_f2].mean(), 1)
+                self.assertLess(delta, 4.0, f"Impulse response estimation for band {band} deviates too much from the real response.")
+
+                # plot the impulse response
+                line = ax.semilogx(ir_res.freq, est_db, label=f"Est.@{(band.center_freq/1000):.1f} kHz (d = {delta} dB)", alpha=0.7)
+                ax.semilogx(freqs, h_db, label=f"Real@{band.center_freq:.1f} Hz", color=line[0].get_color(), linestyle='--', alpha=0.7)
+
+        ax.set_title(f"Impulse Response Estimation (SNR={snr} dB, nfft={nfft}, order={order}, ntaps={ntaps}, use_activity={use_activity})")
+        ax.set_xlabel("Frequency (Hz)")
+        ax.set_ylabel("Magnitude (dB)")
+        ax.legend(ncols=2, bbox_to_anchor=(1.05, 1.0), loc='upper left')
+        ax.grid(True, which='both', linestyle='--', linewidth=0.5)
+        ax.set_ylim(-25, 5)
+        ax.set_xlim(50, 16000)
+        fig.tight_layout()
+        fig.savefig(self.output_dir/f"ir_estimate_snr={int(snr)}db_{nfft}_{order}_{ntaps}_act={use_activity}.svg")
+        #plt.show()
+
+if __name__ == '__main__':
+    unittest.main()
\ No newline at end of file

tests/test_level_calc.py

+import unittest
+import numpy as np
+import soundfile as sf
+from pathlib import Path
+import matplotlib.pyplot as plt
+import pandas as pd
+from typing import Tuple
+from numpy.typing import NDArray
+
+from nonlinearity import MIN_LEVEL_RMS, FS_GLOBAL
+from nonlinearity.spectrum import FftSpectrumVsTime
+from nonlinearity.spectrum.bands import FrequencyBands, get_bands_iec_61260_1_2014
+from nonlinearity.spectrum.weighting import get_weighting
+
+from tests import get_output_dir
+from tests.testfiles.file_dl import get_p501_file
+from tests.testfiles import get_temp_dir
+
+plt_args = {'alpha': 0.7, 'lw': 3.0}
+
+class LevelCalcTestCase(unittest.TestCase):
+    P501_FILE = 'P501_D_AM_fm_FB_48k.wav'  # default file for level tests
+
+    @classmethod
+    def setUpClass(cls):
+        # read the audio file
+        src_file = get_p501_file(cls.P501_FILE, clause='D')
+        if not src_file.is_file():
+            raise FileNotFoundError(f"Test file {src_file} not found.")
+
+        cls._signal: Tuple[NDArray, int] = sf.read(src_file)
+
+    def test_weighting(self, fs: int = FS_GLOBAL, eval_fmin: float = 100.0, eval_fmax: float = 12000.0, fmin: float = 20, fmax: float = 20000):
+        freq = np.fft.rfftfreq(8192, d=1/fs)
+        freq_default = freq[(freq >= fmin) & (freq <= fmax)]
+        ref_dir = get_temp_dir(False, 'ref')
+
+        # check all weightings
+        fig, ax = plt.subplots(figsize=(10, 6))
+        for w in ['A', 'B', 'C', 'D', 'itu_r_468']:
+            with self.subTest(weighting=w):
+                # load ref-data, if available
+                ref_file = ref_dir / f'{w}.txt'
+                if ref_file.is_file():
+                    ref_freq, ref_w_db = np.loadtxt(ref_file, unpack=True)
+                else:
+                    ref_freq = freq_default
+                    ref_w_db = None
+
+                # calculate weighting
+                weighting = get_weighting(ref_freq, w, representation='dB')
+                self.assertIsNotNone(weighting)
+
+                l1 = ax.semilogx(ref_freq, weighting, label=f'{w}-weighting', alpha=0.7)
+                if ref_w_db is not None:
+                    ax.semilogx(ref_freq, ref_w_db, label=f'{w}-weighting (reference)', alpha=0.7, linestyle='--', color=l1[0].get_color())
+                    # frequency range for evaluation
+                    ii_eval = (ref_freq >= eval_fmin) & (ref_freq <= eval_fmax)
+                    self.assertLessEqual(np.abs(ref_w_db-weighting)[ii_eval].max(), 0.30, msg=f"Weighting {w} differs too much from reference in range {eval_fmin}-{eval_fmax} Hz.")
+
+        ax.grid(True)
+        ax.set_title("Weighting References")
+        ax.set_xlabel('Frequency (Hz)')
+        ax.set_ylabel('Magnitude (dB)')
+        ax.legend()
+        fig.tight_layout()
+        fig.savefig(get_output_dir(sub_dir='plots') / 'weighting_vs_references.svg')
+        #plt.show()
+
+
+    def test_bands(self, fs: int = 48000, fmin: float = 20.0, fmax: float = 20000.0, min_order: int = 16):
+        for fraction in [1, 3, 6, 12]:
+            with self.subTest(fraction=fraction):
+                bands = get_bands_iec_61260_1_2014(fraction=fraction, fmin=fmin, fmax=fmax)
+                self.assertIsInstance(bands, FrequencyBands)
+                self.assertGreater(len(bands), 0, msg="No bands created.")
+
+                freq = np.fft.rfftfreq(8192, d=1/fs)
+                resp = np.zeros((len(bands), freq.shape[0]), dtype=float)
+
+                fig, ax = plt.subplots(figsize=(10, 6))
+
+                for i, band in enumerate(bands):
+                    resp[i, :] = band.get_freq_resp(fs, freq, order=max(min_order, 24))
+                    ax.semilogx(freq, 20*np.log10(resp[i, :]), label=f"{band.center_freq} Hz", alpha=0.7)
+
+                # correction for the overall response - should be 0 dB between fmin and fmax
+                overall_db = 20*np.log10(np.maximum(resp.sum(axis=0), 1e-10))  # avoid log(0) issues
+                ax.semilogx(freq, overall_db, 'k--', label="sum", alpha=0.7)
+
+                overall_db_range = overall_db[(freq >= fmin) & (freq <= fmax)]
+                print(overall_db_range.mean(), overall_db_range.std(), overall_db_range.min(), overall_db_range.max())
+                ax.set_title(f"Frequency Response of Bands: {bands.title}")
+                ax.grid(True)
+                ax.legend()
+                ax.set_xlabel('Frequency (Hz)')
+                ax.set_ylabel('Magnitude [dB]')
+                ax.set_ylim(-90, 10)
+                fig.tight_layout()
+                fig.savefig(get_output_dir(sub_dir='plots') / f'1-{fraction}-octave-bands.svg')
+
+        #plt.show()
+
+
+    def test_spectrum_calc(self, fraction: int = 3):
+        src_file = Path(self.P501_FILE)
+        base_name = f"{src_file.stem}_FFT8K_Hann_75"
+        s, fs = self._signal
+
+        # check 3D data / FFT vs time:
+        spec_vs_time = FftSpectrumVsTime.from_signals(s, fs)
+        for all_freqs in [True, False]:
+            data_3d = spec_vs_time.get_values(all_freqs=all_freqs)
+            freq, time = spec_vs_time.get_freqs(all_freqs=all_freqs), spec_vs_time.time
+
+            self.assertEqual(data_3d.shape[0], freq.shape[0])
+            self.assertEqual(data_3d.shape[1], time.shape[0])
+
+        # check 2D data / FFT vs time:
+        spec_avg = spec_vs_time.get_aggregate(aggregate='mean')
+        fft_db = spec_avg.get_values('db')
+        freq = spec_avg.get_freqs()
+        self.assertEqual(fft_db.shape[0], freq.shape[0])
+
+        # convert to bands vs time
+        bands_vs_time = spec_vs_time.to_bands_vs_time(fraction=fraction)
+        self.assertEqual(bands_vs_time.time.shape[0], spec_vs_time.time.shape[0], msg="Time dimension mismatch between bands and FFT vs time.")
+
+        # check 2D data / bands vs time:
+        bands_avg = bands_vs_time.get_aggregate(aggregate='mean')
+        bands_db = bands_avg.get_values('db')
+        bands_freq = bands_avg.get_freqs()
+        self.assertEqual(bands_db.shape[0], bands_freq.shape[0])
+
+        # plot & show reference data of FFT and bands
+        ref_dir = get_temp_dir(False, 'ref')
+        self.assertTrue(ref_dir.is_dir(), msg=f"Reference directory {ref_dir} does not exist.")
+
+        df = pd.DataFrame(columns=['SpectrumType', 'Calculated', 'RefFile', 'Weighting'])
+        df.loc[len(df)] = ['FFT', spec_avg, f"{base_name}", 'Z']
+        df.loc[len(df)] = ['FFT', spec_avg, f"{base_name}_A", 'A']
+        df.loc[len(df)] = ['Bands', bands_avg, f"{base_name}_1-{fraction}Oct", 'Z']
+        df.loc[len(df)] = ['Bands', bands_avg, f"{base_name}_1-{fraction}Oct_A", 'A']
+        #df = df.sort_values(by=['SpectrumType', 'Weighting'], ascending=[True, True])
+        df['RefFile'] = df['RefFile'].apply(lambda x: ref_dir / Path(x).with_suffix('.txt'))
+
+        spec_types = df['SpectrumType'].unique()
+        fig, axes = plt.subplots(spec_types.shape[0], 1, figsize=(10, 6))
+        axes = np.atleast_1d(axes)
+
+        for i, (spec_type, df_st) in enumerate(df.groupby('SpectrumType', sort=False)):
+            ax = axes[i]
+            ax.set_title(f"{spec_type} Spectrum for {src_file.name}")
+            ax.set_xlabel('Frequency (Hz)')
+            ax.set_ylabel('Magnitude (dB)')
+            ax.grid(True)
+
+            for weight, df_weight in df_st.groupby('Weighting', sort=False):
+                for key, row in df_weight.iterrows():
+                    if (ref_file := row['RefFile']).is_file():
+                        ref_freq, ref_spec_db = np.loadtxt(ref_file, unpack=True)
+                        calc_spec_db = row['Calculated'].get_values('db', weight, all_freqs=False)
+                        calc_freq = row['Calculated'].get_freqs(all_freqs=False)
+
+                        l1 = ax.semilogx(calc_freq, calc_spec_db, label=f"Calculated (Weight={weight})", linestyle='-', **plt_args)
+                        ax.semilogx(ref_freq, ref_spec_db, label=f"Reference (Weight={weight})", linestyle='--', color=l1[0].get_color(), **plt_args)
+
+            ax.legend()
+
+        fig.tight_layout()
+        plt.show()
+
+    def test_level_vs_time(self, fmin: float = 100.0, fmax: float = 12000.0):
+        # plot & show reference data of FFT and bands
+        ref_dir = get_temp_dir(False, 'ref')
+        self.assertTrue(ref_dir.is_dir(), msg=f"Reference directory {ref_dir} does not exist.")
+
+        src_file = Path(self.P501_FILE).stem
+        s, fs = self._signal
+        time, ref_lvt = np.loadtxt(ref_dir / f"{src_file}_LvT_SlRect_43ms.txt", unpack=True)
+        time_a, ref_lvt_a = np.loadtxt(ref_dir / f"{src_file}_LvT_SlRect_43ms_A.txt", unpack=True)
+
+        # check 3D data / FFT vs time:
+        spec_vs_time = FftSpectrumVsTime.from_signals(s, fs, fmin=0, fmax=fs/2)
+        lvt = spec_vs_time.get_level_vs_time().get_values('dB', 'Z')
+        lvt_a = spec_vs_time.get_level_vs_time().get_values('dB', 'A')
+        self.assertIsInstance(lvt, np.ndarray, msg="Level vs time should return a numpy array.")
+        self.assertEqual(lvt.shape[0], spec_vs_time.time.shape[0], msg="Time dimension mismatch between level vs time and FFT vs time.")
+
+        # skip first FFT block
+        dt = spec_vs_time.time[0] # - spec_vs_time.time[0]
+        time, time_a = time[time >= dt], time_a[time_a >= dt]  # remove last time point, as it is not valid for level vs time
+        ref_lvt, ref_lvt_a = ref_lvt[(ref_lvt.shape[0] - time.shape[0]):], ref_lvt_a[(ref_lvt_a.shape[0] - time_a.shape[0]):]
+
+        # plot level vs time
+        fig, ax = plt.subplots(figsize=(10, 6))
+        l1 = ax.plot(spec_vs_time.time, lvt, label='Z-weighting', **plt_args)
+        ax.plot(time, ref_lvt, label='Z-weighting (Ref.)', linestyle='--', color=l1[0].get_color(), **plt_args)
+
+        l2 = ax.plot(spec_vs_time.time, lvt_a, label='A-weighting', **plt_args)
+        ax.plot(time_a, ref_lvt_a, label='A-weighting (Ref.)', linestyle='--', color=l2[0].get_color(), **plt_args)
+
+        ax.set_title(f"Level vs Time for {src_file}")
+        ax.set_xlabel('Time (s)')
+        ax.set_ylabel('Level (dB)')
+        ax.grid(True)
+        ax.legend()
+        fig.tight_layout()
+        fig.savefig(get_output_dir(sub_dir='plots') / f'level_vs_time_{src_file}.svg')
+        plt.show()
+
+    def test_level_calc(self, fmin: float = 100.0, fmax: float = 12000.0):
+        #src_file = Path(self.P501_FILE)
+        s, fs = self._signal
+
+        # check 3D data / FFT vs time:
+        spec_vs_time = FftSpectrumVsTime.from_signals(s, fs, fmin=fmin, fmax=fmax)
+        spec_avg = spec_vs_time.get_aggregate(aggregate='mean')
+        bands_avg = spec_vs_time.to_bands_vs_time(fraction=12).get_aggregate(aggregate='mean')
+
+        # check single values
+        ref_levels = {
+            'Z': -27.12,
+            'A': -30.99,
+            'B': -28.16,
+            'C': -27.23,
+        }
+
+        for w, ref_level in ref_levels.items():
+            with self.subTest(weighting=w, mode='spectrum'):
+                level = spec_avg.get_level('db', w, all_freqs=False)
+                self.assertAlmostEqual(np.abs(level - ref_level), 0.0, delta=0.5, msg=f"Level for {w}-weighting differs too much from reference value {ref_level} dB.")
+
+            with self.subTest(weighting=w, mode='bands'):
+                level_bands = bands_avg.get_level('db', w, all_freqs=False)
+                self.assertAlmostEqual(np.abs(level_bands - ref_level), 0.0, delta=0.5, msg=f"Band level for {w}-weighting differs too much from reference value {ref_level} dB.")
+
+if __name__ == '__main__':
+    unittest.main()

tests/test_sine_tone_comparison.py

@@ -10,16 +10,15 @@ from typing import Optional, Any, List, Dict
 from numpy.typing import NDArray
 import matplotlib.pyplot as plt
 import numpy as np
+from itertools import product
 import pandas as pd
-from rootutils import setup_root
 
-setup_root(".")
-
-from nonlinearity.spectrum import calculate_spectrum, calculate_spectrum_vs_time
-from nonlinearity.calculate import calculate_nlm
+from nonlinearity.spectrum import FftSpectrumAvg, FftSpectrumVsTime
+from nonlinearity.calculate import calculate_nlm, NLMResult
 from tools.metrics import thd_plus_n, THDResult
 from tools.harmonics import create_harmonics
 
+from tests import get_output_dir
 from tests.testfiles import get_temp_dir
 
 FS = 48000  # 32000
@@ -31,6 +30,69 @@ TEST_FREQUENCIES: Dict[int, List[float]] = {
 
 class SineToneTestCase(unittest.TestCase):
 
+    @staticmethod
+    def _plot_distortion_vs_freq(res: NLMResult, thd_results: pd.DataFrame, plot_file: Optional[Path] = None):
+        fig, ax = plt.subplots(1, 1, figsize=(8, 4))
+
+        for name, spectrum in res.spectra.asdict().items():
+            lw = 1.0
+            alpha = 1.0
+            if name == "measurement":
+                name = "meas"
+                lw = 3.0
+                alpha = 0.7
+
+            ax.semilogx(res.freq, 20 * np.log10(spectrum.squeeze()), label=name, alpha=alpha, linewidth=lw)
+
+        # nlm vs frequency - looks ugly ... ?
+        #nlm_freq = res.spectra.linear.squeeze() / res.spectra.nonlinear.squeeze()
+        #nlm_freq = convolve(nlm_freq,[.01 for _ in range(100)],mode="same")
+        #ax.semilogx(res.freq, 20 * np.log10(nlm_freq), label='NLM vs Freq.')
+
+        ax.semilogx(thd_results['THD+N'], 'k--', linewidth=3.0, marker='o', alpha=0.5, label='THD+N')
+        ax.semilogx(thd_results['THD'], 'r--', linewidth=3.0, marker='o', alpha=0.5, label='THD')
+        ax.set_xlabel('Frequency [Hz]')
+        ax.set_ylabel('Level [dB]')
+        ax.set_ylim([-105, 70])
+        ax.grid(True)
+        ax.legend()
+        fig.tight_layout()
+        if plot_file is not None:
+            fig.savefig(plot_file)
+        else:
+            plt.show()
+
+        plt.close(fig)
+
+    @staticmethod
+    def _plot_spectrum_vs_time(*signals: NDArray, fs: int = FS, fmin: float = 20.0, fmax: float = 20_000.0,
+                               labels: Optional[List[str]] = None,
+                               plot_file: Optional[Path] = None, **spectrum_params):
+
+        spectra = [FftSpectrumVsTime.from_signals(signal, fs=fs, fmin=fmin, fmax=fmax, **spectrum_params) for signal in signals]
+        labels = labels or [f'Signal {i+1}' for i in range(len(spectra))]
+
+        fig, axes = plt.subplots(1, len(spectra), sharey=True, figsize=(10, 5))
+        for i, (spec, ax) in enumerate(zip(spectra, axes)):
+            data = np.clip(spec.get_values('dB'), a_max=3.0, a_min=-90)
+            ax.pcolormesh(spec.time, spec.get_freqs(), data, label=labels[i], shading='gouraud', rasterized=True)
+            ax.set_yscale('log')
+            ax.set_ylabel('Frequency [Hz]')
+            ax.set_xlabel('Time [s]')
+            ax.legend()
+
+        fig.tight_layout()
+        if plot_file is not None:
+            fig.savefig(plot_file, dpi=250)
+        else:
+            plt.show()
+
+        plt.close(fig)
+
+    @classmethod
+    def setUpClass(cls):
+        cls.plot_dir = get_output_dir('sine_distortion')
+
     @staticmethod
     def _generate_sine(f0: float, duration: float, fs: int = FS, level: float = -3.01,
                        initial_pause: float = 0.0, trailing_pause: float = 0.0,
@@ -63,12 +125,18 @@ class SineToneTestCase(unittest.TestCase):
 
     def _test_multi_sine(self, test_freq_nbr: int, fmin: float, fmax: float, nperseg: int,
                          duration: float, pause, level: float, snr: float):
+        prefix = f"Test{test_freq_nbr:02d}"
 
         harmonics = [2, 3, 4, 5]
-        df = pd.DataFrame(columns=['NLM', 'THD+N', 'THD'])
-        fig, axs = plt.subplots(3, 2, squeeze=True, sharex=True, sharey=True)
-        #print(axs.flatten())
-        for decrease_in_db_per_octave, ax in zip([6.0, 12, 18, 24, 30, 36], axs.flatten()):
+
+        level_decrease_db = np.arange(12, 36+1, 24, dtype=float)  # [6.0, 12, 18, 24, 30, 36]
+        thd_db = np.arange(5, 50, 5, dtype=float)  # [10, 20, 30, 40, 50]
+
+        df = pd.DataFrame(columns=['NLM', 'THD+N', 'THD'] + (idx := ['Target-THD', 'Decrease'])).set_index(idx)
+
+        for (target_thd_db, decrease_in_db_per_octave) in product(thd_db, level_decrease_db):
+            key = target_thd_db, decrease_in_db_per_octave
+
             spectrum_params = dict(nperseg=nperseg, noverlap=round(0.75 * nperseg), window='hann')
             ir = dict(ntaps=nperseg//2, nperseg=nperseg, overlap=0.75)
 
@@ -88,6 +156,7 @@ class SineToneTestCase(unittest.TestCase):
                 h = np.copy(s)
                 h[:round(FS * duration)] = create_harmonics(s[:round(FS * duration)], FS, harmonics=harmonics,
                                                             harmonic_level_offsets_db=level_loss_db,
+                                                            overall_level_harmonics_db=level - target_thd_db,
                                                             harmonics_only=False)
                 sine_tones_dist += [h]
 
@@ -106,7 +175,7 @@ class SineToneTestCase(unittest.TestCase):
             thd_results = pd.DataFrame(columns=['THD+N', 'THD'])
             for i, s in enumerate(slices):
                 result: THDResult = thd_plus_n(multi_sine_deg[s], FS, f0 := TEST_FREQUENCIES[test_freq_nbr][i],
-                                               harmonics=4, fmin=fmin, fmax=fmax,
+                                               harmonics=len(harmonics)+1, fmin=fmin, fmax=fmax,
                                                **spectrum_params)
 
                 thd_results.loc[f0, ['THD+N', 'THD']] = [result.thdn, result.thd]
@@ -115,107 +184,67 @@ class SineToneTestCase(unittest.TestCase):
             spectrum_params['fraction'] = None
             analysis_params = dict(spectrum=spectrum_params, fmin=fmin, fmax=fmax, ir=ir)
 
-            res = calculate_nlm(multi_sine_deg, multi_sine, FS, compensate_delay=False, **analysis_params)
-            df.loc[decrease_in_db_per_octave, ['NLM', 'THD+N', 'THD']] = [res.nlm, thd_results['THD+N'].mean(),
-                                                                          thd_results['THD'].mean()]
+            res = calculate_nlm(multi_sine_deg, multi_sine, fs=FS, compensate_delay=False, **analysis_params)
+            df.loc[key, ['NLM', 'THD+N', 'THD']] = np.array([res.nlm.squeeze(), *thd_results.mean(axis=0)]).round(1)
 
             # plots - distortion vs frequency
-            #fig, ax = plt.subplots(1, 1, figsize=(10, 5))
-
-            for name, spectrum in res.spectra.asdict().items():
-                if name == "measurement": name = "meas"
-                ax.semilogx(res.freq, 20 * np.log10(spectrum.squeeze()), label=name)
-
-            # nlm vs frequency - looks ugly ... ?
-            #nlm_freq = res.spectra.linear.squeeze() / res.spectra.nonlinear.squeeze()
-            #nlm_freq = convolve(nlm_freq,[.01 for _ in range(100)],mode="same")
-            #ax.semilogx(res.freq, 20 * np.log10(nlm_freq), label='NLM vs Freq.')
-
-            ax.semilogx(thd_results['THD+N'], 'k--', linewidth=3.0, marker='o', alpha=0.5, label='THD+N')
-            ax.semilogx(thd_results['THD'], 'r--', linewidth=3.0, marker='o', alpha=0.5, label='THD')
-            ax.set_xlabel('Frequency [Hz]')
-            ax.set_ylabel('Level [dB]')
-            ax.set_ylim([-105, 70])
-            ax.grid(True)
-            ax.legend()
-            fig.tight_layout()
-            fig.savefig(get_temp_dir() / Path(
-                f'Test{test_freq_nbr:02d}_Distortion_decrease={round(decrease_in_db_per_octave):02d}dB.svg'))
+            plot_file = self.plot_dir / Path(f'{prefix}_TargetTHD={round(target_thd_db):02d}dB_Decr={round(decrease_in_db_per_octave):02d}dB.svg')
+            self._plot_distortion_vs_freq(res, thd_results, plot_file=plot_file)
 
-            """"""
             # plots - spectrum vs time
-            spectrum_params['fraction'] = 12
-            freq, time, spec_ref, _ = calculate_spectrum_vs_time(multi_sine, FS,fmin=fmin, fmax=fmax, **spectrum_params)
-            freq, time, spec_deg, _ = calculate_spectrum_vs_time(multi_sine_deg, FS, fmin=fmin, fmax=fmax,
-                                                                  **spectrum_params)
-
-            data_ref = np.clip(20 * np.log10(spec_ref.squeeze()), a_max=3.0, a_min=-90)
-            data_deg = np.clip(20 * np.log10(spec_deg.squeeze()), a_max=3.0, a_min=-90)
-
-            fig, axes = plt.subplots(1, 2, sharey=True, figsize=(10, 5))
-            axes[0].pcolormesh(time, freq, data_ref, label='source', shading='gouraud', rasterized=True)
-            axes[0].set_yscale('log')
-            axes[0].set_ylabel('Frequency [Hz]')
-            axes[0].set_xlabel('Time [s]')
+            #spectrum_params['fraction'] = 12
+            #plot_file = self.plot_dir / Path(f'{prefix}_Ref_vs_Deg_decrease={round(decrease_in_db_per_octave):02d}dB.svg')
+            #self._plot_spectrum_vs_time(multi_sine, multi_sine_deg, fs=FS, labels=['Source', 'Clipped'], plot_file=plot_file,
+            #                            fmin=fmin, fmax=fmax, **spectrum_params)
 
-            axes[1].pcolormesh(time, freq, data_deg, label='clipped', shading='gouraud', rasterized=True)
-            axes[1].set_xlabel('Time [s]')
-
-            fig.tight_layout()
-            fig.savefig(get_temp_dir() / Path(f'Ref_vs_Deg_decrease={round(decrease_in_db_per_octave):02d}dB.svg'),
-                        dpi=250, )
-
-        fig.set_size_inches(8, 10)
-        fig.subplots_adjust(bottom=0.075, top=0.97, right=0.97)
-        #plt.show()
 
         plt.close('all')
-        print(df)
-        df.to_excel(get_temp_dir() / f"Test{test_freq_nbr:02d}_ResultValues.xlsx")
+        df.to_excel(self.plot_dir / f"{prefix}_ResultValues.xlsx")
 
     def test_generate_harmonics(self):
         # test with sine tones at different frequencies
         f0 = [80, 315, 1000.0]
         duration = 3.0  # per sine tone
-        level_db = -10.0
+        level_db = +10.0
         snr = 50.0
         trailing_pause = 0.5
 
         # harmonic "configuration":
-        num_harmonics = 4
+        num_harmonics = len(harmonics := [3, 4, 5])
         decrease_in_db_per_octave = 18.0
 
         # spectral analysis
-        spec_params = dict(nperseg=(nperseg := 8192 * 2), noverlap=round(0.75 * nperseg), window='hann')
+        spec_params = dict(nperseg=(nperseg := 2*8192), nfft=nperseg, noverlap=round(0.75 * nperseg), window='hamming')
 
         # generate sine tones and harmonics
         sine = []
-        harmonics = []
+        harmonic_signals = []
         for i in range(n_channels := len(f0)):
             # create w/o any noise and pause for generation of harmonics
-            s = self._generate_sine(f0[i], duration, level=level_db, snr=None, trailing_pause=0.0)
-
-            harmonic_freqs = [f0[i] * (j + 2) for j in range(num_harmonics)]
+            harmonic_freqs = [f0[i] * h_order for h_order in harmonics]
             level_loss_db = [-((fh / f0[i]) - 1) * decrease_in_db_per_octave for fh in harmonic_freqs]
-            h = create_harmonics(s, FS, harmonics=num_harmonics, harmonic_level_offsets_db=level_loss_db)
-            harmonics.append(h)
 
+            s = self._generate_sine(f0[i], duration, level=level_db, snr=None, trailing_pause=0.0)
+            h = create_harmonics(s, FS, harmonics=harmonics, harmonic_level_offsets_db=level_loss_db, harmonics_only=True)
+            harmonic_signals.append(h)
+
+            # reference sine tone with noise and trailing pause:
             sine.append(self._generate_sine(f0[i], duration, level=level_db, snr=snr, trailing_pause=trailing_pause))
 
         sine = np.vstack(sine).T
-        harmonics = np.vstack(harmonics).T
+        harmonic_signals = np.vstack(harmonic_signals).T
 
         # spectral analysis
-        freq, spec_avg_harmonics, noise_est_output1, wcf = calculate_spectrum(harmonics, FS, **spec_params)
-        _, spec_avg_sines, noise_est_output2, wcf = calculate_spectrum(sine, FS, **spec_params)
+        spec_avg_sines = FftSpectrumAvg.from_signals(sine, fs=FS, **spec_params)
+        spec_avg_harmonics = FftSpectrumAvg.from_signals(harmonic_signals, fs=FS, **spec_params)
+        freq = spec_avg_sines.get_freqs()
 
         # plot results
         fig, axes = plt.subplots(n_channels, 1, figsize=(n_channels * 4, 6), sharex=True)
         for i in range(n_channels):
-            expected_harmonics = [(h + 2) * f0[i] for h in range(num_harmonics)]
-            axes[i].semilogx(freq, 20 * np.log10(spec_avg_sines[i]) + wcf, label=f'sine tone f0={f0[i]}')
-            axes[i].semilogx(freq, hspec := 20 * np.log10(spec_avg_harmonics[i]) + wcf,
-                             label=f'harmonics {expected_harmonics}')
+            expected_harmonics = [h_order * f0[i] for h_order in harmonics]
+            axes[i].semilogx(freq, spec_avg_sines.get_values('dB', include_window_correction=True)[:, i], label=f'sine tone f0={f0[i]}')
+            axes[i].semilogx(freq, hspec := spec_avg_harmonics.get_values('dB', include_window_correction=True)[:, i], label=f'harmonics {expected_harmonics}')
 
             for exh in expected_harmonics:
                 f_idx = np.argmin(np.abs(freq - exh))

tests/test_swept_sine_comparison.py

@@ -5,87 +5,70 @@ Created on Sep 19 2024 09:39
 @author: Jan.Reimes
 """
 import unittest
-from typing import List
+from pathlib import Path
+from typing import List, Dict, Any, Optional
 import matplotlib.pyplot as plt
 from matplotlib.axes import Axes
 import numpy as np
 from numpy.typing import NDArray
 from scipy.signal import freqz, lfilter
+from ddt import ddt, data
 
-from tools.farina.sweep import generate_farina_sweep
-from tools.farina.deconvolve import estimate_imp_resp
-from nonlinearity.calculate import calculate_nlm
-from nonlinearity.spectrum import calculate_spectrum, calculate_spectrum_vs_time
+from nonlinearity.calculate import calculate_nlm, NLMResult
+from nonlinearity.spectrum.bands import fft_to_octave_bands, get_bands_iec_61260_1_2014
+from nonlinearity.spectrum import ISpectrumVsTime, FftSpectrumAvg, FftSpectrumVsTime, to_db
 
+from tools.plot_spectra import NLMPlot
+from tools.farina.sweep import generate_farina_sweep, FarinaSweep
+from tools.farina.deconvolve import estimate_imp_resp, DeconvolutionResult
+
+from tests import get_output_dir
 from tests.test_base import BaseTestCase, get_config
 
 AxesArray = NDArray[Axes]
 
+def plot_spectrum_vs_time(spec_vs_time: Dict[str, ISpectrumVsTime], nbr_channels: int, nlm: Optional[NLMResult] = None,
+                          plot_file: Optional[Path] = None):
 
-class SweepsTestCase(BaseTestCase):
-
-    def test_sweep_distortion(self, duration: float = 20, nbr_channels: int = 1, trailing_silence: float = 1.0,
-                              level_db: float = -3.01, delay: float = 0.04, harmonics: List[int] = [2, 3, 4]):
-
-        with get_config(return_hydra_config=False) as cfg:
-            fs = cfg.fs
-            max_nbr_ir_samples = cfg.ir.ntaps
-
-            cfg = cfg.copy() # OmegaConf.to_container(cfg, resolve=True, )
-            cfg['file_measured'] = ''
-            cfg['file_ref'] = ''
-            cfg.spectrum.fraction = 12
-            
-        # source sweep:
-        sweep = generate_farina_sweep(fs, duration, level_db=level_db, fmin=35.0, fmax=16_000,
-                                      trailing_silence=trailing_silence)
-
-        # add harmonic distortions:
-        distorted, harmonic_distortions = self._distort_signal(sweep.signal, sweep.fs, delay, harmonics=harmonics,
-                                                               nbr_channels=nbr_channels,
-                                                               harmonic_level_offsets_db=-6,
-                                                               overall_level_harmonics_db=-20.0,
-                                                               min_harmonic_freq=200.0, max_harmonic_freq=6000.0,
-                                                               )
-        sweep_signal, distorted, harmonic_distortions = self._pad_same_length(sweep.signal, distorted, harmonic_distortions)
-        lin_distorted = distorted - harmonic_distortions
-
-        # calculate nlm
-        nlm = calculate_nlm(distorted, sweep_signal, **cfg)
-
-        # apply deconvolution - calculate one more order than the highest harmonic (which should then be approx. zero)
-        ir_res = estimate_imp_resp(distorted, sweep, max_nbr_ir_samples=max_nbr_ir_samples,
-                                   max_nl_order=max(harmonics) + 2,
-                                   fade_in=True, fade_out=True, fade_window='flattop')
-
-        # spectra of sweep signals
-        spec_vs_time = {'source': sweep_signal, 'distorted': distorted, 'harmonic_distortions': harmonic_distortions, 'lin_distorted': lin_distorted}
-        levels = dict()
-        for key, s in spec_vs_time.items():
-            freq, time, spec, wcf = calculate_spectrum_vs_time(s, fs, **cfg.spectrum)
-            spec_avg = np.abs(spec).mean(axis=1)
-            spec_avg[(freq > 0) & (freq < fs/2)] *= np.sqrt(2)
-            levels[key] = np.around(10 * np.log10(np.sum(spec_avg ** 2, axis=0)) + wcf, 2)
-
-            spec_vs_time[key] = self.fr_to_db(spec.squeeze())
-
-        # plot sweep
     fig, axes = plt.subplots(n_plots := len(spec_vs_time), nbr_channels, sharey=True, figsize=(12, 2*n_plots))
     axes: AxesArray = np.atleast_2d(axes).reshape(n_plots, nbr_channels)
     for j, (key, spec) in enumerate(spec_vs_time.items()):
         axes[j, 0].set_ylabel('Frequency [Hz]')
+        time, freq = spec.time, spec.get_freqs()
+        sp_data = spec.get_values('dB')
+        levels = np.around(np.atleast_1d(spec.get_aggregate().get_level('dB')))
+
         for i in range(nbr_channels):
             ax = axes[j, i]
-                ax.pcolormesh(time, freq, spec if nbr_channels == 1 else spec[:,:,i], label=f"{key} Ch{i+1}", shading='gouraud', rasterized=True)
-                ax.set_title(f"{key} {levels[key][i]} dB")
+
+            ax.pcolormesh(time, freq, sp_data[:, :, i], label=f"{key} Ch{i+1}", shading='gouraud', rasterized=True)
+            ax.set_title(f"{key} {levels[i]} dB")
             ax.set_yscale('log')
-                ax.set_ylim([cfg.fmin, cfg.fmax])
+            #ax.set_ylim([cfg.fmin, cfg.fmax])
 
     [ax.set_xlabel('Time [s]') for ax in axes[-1, :]]
-        fig.suptitle(f"NLM = {nlm.nlm[0]:.1f} dB")
+    if isinstance(nlm, NLMResult):
+        for i in range(nbr_channels):
+            title = axes[0, i].get_title()
+            title = f"{title} - NLM = {nlm.nlm[i]:.1f} dB"
+            axes[0, i].set_title(title)
+
     fig.tight_layout()
+    if plot_file:
+        fig.savefig(plot_file, bbox_inches='tight')
+    else:
         plt.show()
-        """ """
+
+    plt.close(fig)
+
+def plot_nonlinear_deconvolution(sweep: FarinaSweep, ir_res: DeconvolutionResult, nbr_channels: int,
+                                 nlm: Optional[NLMResult] = None, plot_file: Optional[Path] = None,
+                                 **spec_args: Dict[str, Any]
+                                 ):
+    fs = ir_res.fs
+    fraction = spec_args.get('fraction')
+    bands = get_bands_iec_61260_1_2014(fraction=spec_args.get('fraction', 24), fmin=spec_args.get('fmin', 20.0),
+                                        fmax=spec_args.get('fmax', 20000.0))
 
     # plot linear & non-linear impulse responses/frequency responses from Farina deconvolution
     fig, axes = plt.subplots(2, nbr_channels)
@@ -96,12 +79,11 @@ class SweepsTestCase(BaseTestCase):
 
         # frequency responses of linear and non-linear components
         for nl_order, imp_resp in ir_res.nl_responses.items():
-                freq, fr = freqz(ir := imp_resp[:, i], fs=fs, worN=cfg.spectrum.nperseg)
+            freq, fr = freqz(ir := imp_resp[:, i], fs=fs, worN=spec_args.get('nperseg', 1024))
 
-                if cfg.spectrum.fraction:
-                    fr, freq = self.fft_to_octave_bands(fr, fs, fraction=cfg.spectrum.fraction, is_transfer_function=True)
+            if fraction:
+                freq, fr = fft_to_octave_bands(fr, fs, bands, input_domain='rms', output_domain='dB', is_transfer_function=True)
 
-                fr = self.fr_to_db(fr)
             axes[0, i].semilogx(freq, fr, label=f"linear Ch{i + 1}" if nl_order < 2 else f"harmonic#{nl_order} Ch{i + 1}")
 
     # estimate signal/spectrum of NL components
@@ -113,27 +95,105 @@ class SweepsTestCase(BaseTestCase):
 
             nl_signal_est[:, i] += lfilter(imp_resp[:, i], [1.0], sweep.signal)
 
-            freq, nl_spec_est, _, _ = calculate_spectrum(nl_signal_est[:, i], fs=fs, **cfg.spectrum)
-            axes[1, i].semilogx(freq, self.fr_to_db(nl_spec_est.squeeze()), label=f"Sum(Harmonics) Ch{i + 1}", alpha=0.5)
-
+        #freq, nl_spec_est, _, _ = calculate_spectrum(nl_signal_est[:, i], fs=fs, **cfg.spectrum)
+        nl_spec_est = FftSpectrumAvg.from_signals(nl_signal_est[:, i], fs, **spec_args)
 
-            nl_pow2 = 10 * np.log10()
+        axes[1, i].semilogx(nl_spec_est.get_freqs(), nl_spec_est.get_values('dB'), label=f"Sum(Harmonics) Ch{i + 1}", alpha=0.5)
 
     # from NLM
+    """
     for i, ir in enumerate(nlm.freq_resp_result.imp_resp.T):
-            fr_lin = self.fr_to_db(nlm.freq_resp_result.freq_resp[:, i])
+        fr_lin = to_db(nlm.freq_resp_result.freq_resp[:, i], True)
         axes[0, i].semilogx(nlm.freq_resp_result.freq, fr_lin, label=f"IR(NLM) Ch{i + 1}", alpha=0.5, linewidth=4.0)
 
-            nl_spec = self.fr_to_db(nlm.spectra.nonlinear[:, i])
+        nl_spec = to_db(nlm.spectra.nonlinear[:, i], True)
         axes[1, i].semilogx(nlm.freq, nl_spec, label=f"NLM Ch{i + 1}", alpha=0.5, linewidth=4.0)
+    """
 
     for ax in axes.flatten():
         ax.grid(True)
         ax.legend()
 
     fig.tight_layout()
+    if plot_file:
+        fig.savefig(plot_file, bbox_inches='tight')
+    else:
         plt.show()
 
+    plt.close(fig)
+
+
+@ddt
+class SweepsTestCase(BaseTestCase):
+
+    @classmethod
+    def setUpClass(cls):
+        super().setUpClass()
+        cls.plot_dir = get_output_dir('plots_sweeps')
+
+    @data(-40, -30, -20, -10)
+    def test_sweep_distortion(self, level_harmonics_db: float, duration: float = 20, nbr_channels: int = 1, trailing_silence: float = 1.0,
+                              src_level_db: float = -3.01, delay: float = 0.04, harmonics: List[int] = [2, 3]):
+
+        with get_config(return_hydra_config=False) as cfg:
+            fs = cfg.fs #= 48000
+            max_nbr_ir_samples = cfg.ir.ntaps
+
+            cfg = cfg.copy() # OmegaConf.to_container(cfg, resolve=True, )
+            cfg['file_measured'] = ''
+            cfg['file_ref'] = ''
+            cfg.ir.nperseg = 16384
+            cfg.ir.window = 'hann'
+            #cfg.spectrum.fraction = 24
+
+
+        # source sweep:
+        sweep: FarinaSweep = generate_farina_sweep(fs, duration, level_db=src_level_db, fmin=35.0, fmax=16_000,
+                                      trailing_silence=trailing_silence)
+
+        # add harmonic distortions:
+        distorted, harmonic_distortions = self._distort_signal(sweep.signal, sweep.fs, delay, harmonics=harmonics,
+                                                               nbr_channels=nbr_channels,
+                                                               harmonic_level_offsets_db=-6*np.arange(1, len(harmonics)+1),
+                                                               overall_level_harmonics_db=level_harmonics_db,
+                                                               min_harmonic_freq=200.0, max_harmonic_freq=6000.0,
+                                                               #rir_filter_id=None
+                                                               )
+        sweep_signal, distorted, harmonic_distortions = self._pad_same_length(sweep.signal, distorted, harmonic_distortions)
+        lin_distorted = distorted - harmonic_distortions
+
+        # calculate nlm
+        nlm = calculate_nlm(distorted, sweep_signal, **cfg)
+        #NLMPlot(nlm, title='NLM', params=cfg).plot(plot_components=True, plot_nlm=True, plot_fr=True, plot_cxy=True, show_blocking=True)
+
+        # apply deconvolution - calculate one more order than the highest harmonic (which should then be approx. zero)
+        ir_res = estimate_imp_resp(distorted, sweep, max_nbr_ir_samples=max_nbr_ir_samples,
+                                   max_nl_order=max(harmonics) + 1,
+                                   fade_in=True, fade_out=True, # fade_window='flattop'
+                                   )
+
+        # spectra of sweep signals
+        spec_vs_time = {'source': sweep_signal, 'distorted': distorted, 'harmonic_distortions': harmonic_distortions, 'lin_distorted': lin_distorted}
+        levels = dict()
+        for key, s in spec_vs_time.items():
+            spec = FftSpectrumVsTime.from_signals(s, fs, always_2d=True, **cfg.spectrum)
+            levels[key] = [np.around(spec.get_aggregate().get_level('dB'), 2)]
+            spec_vs_time[key] = spec
+
+            #freq, time, spec, wcf = calculate_spectrum_vs_time(s, fs, **cfg.spectrum)
+            #levels[key] = [np.around(calculate_level(s, fs, **cfg.spectrum), 2)]
+            #spec_vs_time[key] = self.fr_to_db(spec.squeeze())
+
+        # plot deconvolution
+        plot_nonlinear_deconvolution(sweep, ir_res, nbr_channels,
+                                     nlm=nlm,
+                                     plot_file=self.plot_dir / f'sweep_nonlinearities_{level_harmonics_db}dB.svg',
+                                     **cfg.spectrum
+                                     )
+
+        # plot sweep
+        plot_spectrum_vs_time(spec_vs_time, nbr_channels, nlm=nlm, plot_file=self.plot_dir / f'sweep_distortion_{level_harmonics_db}dB.svg')
+
 
 if __name__ == '__main__':
     unittest.main()

tests/testfiles/file_dl.py

@@ -6,13 +6,15 @@ Created on Jun 29 2024 00:11
 """
 
 from pathlib import Path
+from typing import LiteralString
 
 from tools.downloader import Downloads, FileFetcher
 from tests.testfiles import get_download_dir, get_extract_dir
 
 
-def get_p501_file(name: str) -> Path:
-    return FileFetcher(Downloads.P501, dl_dir=get_download_dir(), extract_dir=get_extract_dir()).get_file(name)
+def get_p501_file(name: str, clause: LiteralString|None = None) -> Path:
+    dl = Downloads.P501 if clause is None else Downloads.get(f'P501_{clause.upper()}')
+    return FileFetcher(dl, dl_dir=get_download_dir(), extract_dir=get_extract_dir()).get_file(name)
 
 
 def get_g191_file(name: str) -> Path:

tests/testfiles/ref/P501_D_AM_fm_FB_48k_FFT8K_Hann_75.txt

+5.859375	-111.86
+11.71875	-105.25
+17.578125	-99.06
+23.4375	-89.96
+29.296875	-82.39
+35.15625	-78.02
+41.015625	-74.25
+46.875	-71.70
+52.734375	-68.77
+58.59375	-66.12
+64.453125	-64.21
+70.3125	-62.53
+76.171875	-62.94
+82.03125	-64.27
+87.890625	-64.27
+93.75	-56.13
+99.609375	-51.07
+105.46875	-49.15
+111.328125	-43.33
+117.1875	-40.32
+123.046875	-42.40
+128.90625	-43.91
+134.765625	-43.02
+140.625	-43.65
+146.484375	-45.84
+152.34375	-43.90
+158.203125	-41.88
+164.0625	-41.38
+169.921875	-41.24
+175.78125	-43.12
+181.640625	-46.26
+187.5	-47.28
+193.359375	-48.80
+199.21875	-49.87
+205.078125	-52.11
+210.9375	-53.19
+216.796875	-53.23
+222.65625	-49.94
+228.515625	-44.71
+234.375	-45.05
+240.234375	-46.16
+246.09375	-46.90
+251.953125	-50.70
+257.8125	-52.64
+263.671875	-51.53
+269.53125	-49.63
+275.390625	-48.06
+281.25	-48.58
+287.109375	-50.83
+292.968751	-51.06
+298.828126	-50.27
+304.687501	-49.12
+310.546876	-48.87
+316.406251	-46.44
+322.265626	-44.72
+328.125001	-45.02
+333.984376	-45.13
+339.843751	-44.58
+345.703126	-44.73
+351.562501	-46.63
+357.421876	-47.92
+363.281251	-47.51
+369.140626	-48.36
+375.000001	-50.00
+380.859376	-51.13
+386.718751	-52.37
+392.578126	-52.52
+398.437501	-51.68
+404.296876	-51.75
+410.156251	-50.69
+416.015626	-50.08
+421.875001	-51.56
+427.734376	-53.33
+433.593751	-52.17
+439.453126	-52.36
+445.312501	-52.56
+451.171876	-49.47
+457.031251	-44.62
+462.890626	-44.93
+468.750001	-47.04
+474.609376	-44.62
+480.468751	-41.35
+486.328126	-39.07
+492.187501	-40.48
+498.046876	-43.56
+503.906251	-45.56
+509.765626	-48.19
+515.625001	-49.23
+521.484376	-48.99
+527.343751	-49.60
+533.203126	-50.51
+539.062501	-49.85
+544.921876	-48.13
+550.781251	-45.50
+556.640626	-44.38
+562.500001	-44.90
+568.359376	-46.07
+574.218751	-45.50
+580.078126	-45.66
+585.937501	-46.09
+591.796876	-46.23
+597.656251	-46.50
+603.515626	-46.77
+609.375001	-46.55
+615.234376	-47.01
+621.093751	-47.54
+626.953126	-48.82
+632.812501	-49.77
+638.671876	-50.37
+644.531251	-51.24
+650.390626	-52.41
+656.250001	-52.42
+662.109376	-52.55
+667.968751	-51.22
+673.828126	-49.00
+679.687501	-47.62
+685.546876	-45.79
+691.406251	-44.34
+697.265626	-45.19
+703.125001	-45.75
+708.984376	-45.92
+714.843751	-46.84
+720.703126	-48.29
+726.562501	-48.63
+732.421876	-48.90
+738.281251	-48.94
+744.140626	-49.67
+750.000001	-51.22
+755.859376	-52.26
+761.718751	-54.04
+767.578126	-55.41
+773.437501	-56.39
+779.296876	-58.02
+785.156251	-59.01
+791.015626	-57.88
+796.875001	-56.90
+802.734376	-56.45
+808.593751	-56.03
+814.453126	-56.40
+820.312501	-54.98
+826.171876	-53.84
+832.031251	-54.27
+837.890626	-54.63
+843.750001	-52.33
+849.609376	-52.02
+855.468751	-52.30
+861.328126	-51.31
+867.187501	-52.11
+873.046876	-53.40
+878.906252	-55.73
+884.765627	-57.52
+890.625002	-59.49
+896.484377	-59.74
+902.343752	-60.35
+908.203127	-60.09
+914.062502	-59.84
+919.921877	-59.43
+925.781252	-57.62
+931.640627	-58.05
+937.500002	-57.53
+943.359377	-57.76
+949.218752	-59.31
+955.078127	-59.57
+960.937502	-57.58
+966.796877	-56.18
+972.656252	-56.55
+978.515627	-59.35
+984.375002	-59.17
+990.234377	-57.71
+996.093752	-58.43
+1001.953127	-59.99
+1007.812502	-59.30
+1013.671877	-59.45
+1019.531252	-58.97
+1025.390627	-59.09
+1031.250002	-58.66
+1037.109377	-57.23
+1042.968752	-57.02
+1048.828127	-58.35
+1054.687502	-59.45
+1060.546877	-60.66
+1066.406252	-61.89
+1072.265627	-60.65
+1078.125002	-61.55
+1083.984377	-61.55
+1089.843752	-60.01
+1095.703127	-57.75
+1101.562502	-55.82
+1107.421877	-56.36
+1113.281252	-57.75
+1119.140627	-59.42
+1125.000002	-56.38
+1130.859377	-57.63
+1136.718752	-58.71
+1142.578127	-60.89
+1148.437502	-60.10
+1154.296877	-59.06
+1160.156252	-57.32
+1166.015627	-57.25
+1171.875002	-56.94
+1177.734377	-56.84
+1183.593752	-58.98
+1189.453127	-60.97
+1195.312502	-61.18
+1201.171877	-61.34
+1207.031252	-60.80
+1212.890627	-61.62
+1218.750002	-62.65
+1224.609377	-63.37
+1230.468752	-64.37
+1236.328127	-65.23
+1242.187502	-66.00
+1248.046877	-66.36
+1253.906252	-66.29
+1259.765627	-64.60
+1265.625002	-63.29
+1271.484377	-60.90
+1277.343752	-59.59
+1283.203127	-58.56
+1289.062502	-57.97
+1294.921877	-58.79
+1300.781252	-59.96
+1306.640627	-60.31
+1312.500002	-61.48
+1318.359377	-63.95
+1324.218752	-64.89
+1330.078127	-63.96
+1335.937502	-61.78
+1341.796877	-61.47
+1347.656252	-62.25
+1353.515627	-61.89
+1359.375002	-61.52
+1365.234377	-63.82
+1371.093752	-65.47
+1376.953127	-64.39
+1382.812502	-63.47
+1388.671877	-61.28
+1394.531252	-60.64
+1400.390627	-61.45
+1406.250002	-61.53
+1412.109377	-63.20
+1417.968752	-65.53
+1423.828127	-68.49
+1429.687502	-70.08
+1435.546877	-69.18
+1441.406252	-68.98
+1447.265627	-67.80
+1453.125002	-67.10
+1458.984377	-66.98
+1464.843752	-66.88
+1470.703128	-67.44
+1476.562503	-68.26
+1482.421878	-67.97
+1488.281253	-68.29
+1494.140628	-67.83
+1500.000003	-64.56
+1505.859378	-62.66
+1511.718753	-62.66
+1517.578128	-63.76
+1523.437503	-64.45
+1529.296878	-66.21
+1535.156253	-67.55
+1541.015628	-68.33
+1546.875003	-68.46
+1552.734378	-69.43
+1558.593753	-68.52
+1564.453128	-67.36
+1570.312503	-65.98
+1576.171878	-64.70
+1582.031253	-64.24
+1587.890628	-64.50
+1593.750003	-65.06
+1599.609378	-65.31
+1605.468753	-64.85
+1611.328128	-64.66
+1617.187503	-63.22
+1623.046878	-64.24
+1628.906253	-66.86
+1634.765628	-66.76
+1640.625003	-66.68
+1646.484378	-66.36
+1652.343753	-65.24
+1658.203128	-65.29
+1664.062503	-64.90
+1669.921878	-63.39
+1675.781253	-63.72
+1681.640628	-64.63
+1687.500003	-64.47
+1693.359378	-64.07
+1699.218753	-62.98
+1705.078128	-60.50
+1710.937503	-59.94
+1716.796878	-61.80
+1722.656253	-64.01
+1728.515628	-63.92
+1734.375003	-64.53
+1740.234378	-65.27
+1746.093753	-65.63
+1751.953128	-66.42
+1757.812503	-66.40
+1763.671878	-66.12
+1769.531253	-64.97
+1775.390628	-65.46
+1781.250003	-64.98
+1787.109378	-64.75
+1792.968753	-64.81
+1798.828128	-64.80
+1804.687503	-65.32
+1810.546878	-66.30
+1816.406253	-64.33
+1822.265628	-63.35
+1828.125003	-64.13
+1833.984378	-65.48
+1839.843753	-64.54
+1845.703128	-64.72
+1851.562503	-67.38
+1857.421878	-66.62
+1863.281253	-64.41
+1869.140628	-65.78
+1875.000003	-68.69
+1880.859378	-64.37
+1886.718753	-62.48
+1892.578128	-64.60
+1898.437503	-65.24
+1904.296878	-64.06
+1910.156253	-65.74
+1916.015628	-66.64
+1921.875003	-66.84
+1927.734378	-69.57
+1933.593753	-68.94
+1939.453128	-67.90
+1945.312503	-66.94
+1951.171878	-66.53
+1957.031253	-67.72
+1962.890628	-68.67
+1968.750003	-67.81
+1974.609378	-67.00
+1980.468753	-67.84
+1986.328128	-67.95
+1992.187503	-66.90
+1998.046878	-67.16
+2003.906253	-69.58
+2009.765628	-69.48
+2015.625003	-67.35
+2021.484378	-67.78
+2027.343753	-68.51
+2033.203128	-68.88
+2039.062503	-68.67
+2044.921878	-68.61
+2050.781253	-68.76
+2056.640629	-68.64
+2062.500004	-68.58
+2068.359379	-67.43
+2074.218754	-69.50
+2080.078129	-68.86
+2085.937504	-67.81
+2091.796879	-68.99
+2097.656254	-66.15
+2103.515629	-65.80
+2109.375004	-68.87
+2115.234379	-69.43
+2121.093754	-67.36
+2126.953129	-68.66
+2132.812504	-69.95
+2138.671879	-69.17
+2144.531254	-70.37
+2150.390629	-71.79
+2156.250004	-72.47
+2162.109379	-71.83
+2167.968754	-71.64
+2173.828129	-72.13
+2179.687504	-72.55
+2185.546879	-71.17
+2191.406254	-69.55
+2197.265629	-69.24
+2203.125004	-68.32
+2208.984379	-69.28
+2214.843754	-68.75
+2220.703129	-67.90
+2226.562504	-68.81
+2232.421879	-68.78
+2238.281254	-69.53
+2244.140629	-69.85
+2250.000004	-70.05
+2255.859379	-69.87
+2261.718754	-70.47
+2267.578129	-71.51
+2273.437504	-71.36
+2279.296879	-72.20
+2285.156254	-72.01
+2291.015629	-71.47
+2296.875004	-71.27
+2302.734379	-71.19
+2308.593754	-71.81
+2314.453129	-72.83
+2320.312504	-72.28
+2326.171879	-70.48
+2332.031254	-71.03
+2337.890629	-70.18
+2343.750004	-68.22
+2349.609379	-68.35
+2355.468754	-68.56
+2361.328129	-69.66
+2367.187504	-68.89
+2373.046879	-69.03
+2378.906254	-70.15
+2384.765629	-72.08
+2390.625004	-72.39
+2396.484379	-71.45
+2402.343754	-72.21
+2408.203129	-73.78
+2414.062504	-72.56
+2419.921879	-72.58
+2425.781254	-72.50
+2431.640629	-71.45
+2437.500004	-70.80
+2443.359379	-70.02
+2449.218754	-71.24
+2455.078129	-72.27
+2460.937504	-71.85
+2466.796879	-72.37
+2472.656254	-74.80
+2478.515629	-75.10
+2484.375004	-73.47
+2490.234379	-74.02
+2496.093754	-73.05
+2501.953129	-74.35
+2507.812504	-74.87
+2513.671879	-73.97
+2519.531254	-72.01
+2525.390629	-71.69
+2531.250004	-71.84
+2537.109379	-72.50
+2542.968754	-71.11
+2548.828129	-71.47
+2554.687504	-72.43
+2560.546879	-70.14
+2566.406254	-69.98
+2572.265629	-71.42
+2578.125004	-71.57
+2583.984379	-72.51
+2589.843754	-72.39
+2595.703129	-72.46
+2601.562504	-73.18
+2607.421879	-72.37
+2613.281254	-72.51
+2619.140629	-72.24
+2625.000004	-74.12
+2630.859379	-74.30
+2636.718754	-73.62
+2642.57813	-75.06
+2648.437505	-75.06
+2654.29688	-73.58
+2660.156255	-73.64
+2666.01563	-73.88
+2671.875005	-74.72
+2677.73438	-74.21
+2683.593755	-73.76
+2689.45313	-72.92
+2695.312505	-72.53
+2701.17188	-74.35
+2707.031255	-75.86
+2712.89063	-74.64
+2718.750005	-73.74
+2724.60938	-75.01
+2730.468755	-74.84
+2736.32813	-73.71
+2742.187505	-73.72
+2748.04688	-72.36
+2753.906255	-72.05
+2759.76563	-72.39
+2765.625005	-72.39
+2771.48438	-72.35
+2777.343755	-74.36
+2783.20313	-75.75
+2789.062505	-76.19
+2794.92188	-76.19
+2800.781255	-75.87
+2806.64063	-76.51
+2812.500005	-75.92
+2818.35938	-74.64
+2824.218755	-73.14
+2830.07813	-72.86
+2835.937505	-75.04
+2841.79688	-73.47
+2847.656255	-71.76
+2853.51563	-72.56
+2859.375005	-73.38
+2865.23438	-74.88
+2871.093755	-73.04
+2876.95313	-71.80
+2882.812505	-72.40
+2888.67188	-74.87
+2894.531255	-75.72
+2900.39063	-75.65
+2906.250005	-75.95
+2912.10938	-75.03
+2917.968755	-72.82
+2923.82813	-71.31
+2929.687505	-71.31
+2935.54688	-72.36
+2941.406255	-73.84
+2947.26563	-75.66
+2953.125005	-76.44
+2958.98438	-76.37
+2964.843755	-76.92
+2970.70313	-75.30
+2976.562505	-74.30
+2982.42188	-75.42
+2988.281255	-74.32
+2994.14063	-73.84
+3000.000005	-72.38
+3005.85938	-71.02
+3011.718755	-71.59
+3017.57813	-73.45
+3023.437505	-75.20
+3029.29688	-75.44
+3035.156255	-74.09
+3041.01563	-74.36
+3046.875005	-74.32
+3052.73438	-72.72
+3058.593755	-73.39
+3064.45313	-73.43
+3070.312505	-73.43
+3076.17188	-73.92
+3082.031255	-73.77
+3087.89063	-73.63
+3093.750005	-74.00
+3099.60938	-75.13
+3105.468755	-75.50
+3111.32813	-73.42
+3117.187505	-74.07
+3123.04688	-76.96
+3128.906255	-77.59
+3134.76563	-77.03
+3140.625005	-76.99
+3146.48438	-75.79
+3152.343755	-75.51
+3158.20313	-75.97
+3164.062505	-73.47
+3169.92188	-73.43
+3175.781255	-73.73
+3181.64063	-73.81
+3187.500005	-74.96
+3193.35938	-74.97
+3199.218755	-73.55
+3205.07813	-72.84
+3210.937505	-72.67
+3216.79688	-72.35
+3222.656256	-73.24
+3228.515631	-72.45
+3234.375006	-70.66
+3240.234381	-71.28
+3246.093756	-72.72
+3251.953131	-73.34
+3257.812506	-74.06
+3263.671881	-74.21
+3269.531256	-73.41
+3275.390631	-73.11
+3281.250006	-73.14
+3287.109381	-73.42
+3292.968756	-73.17
+3298.828131	-71.58
+3304.687506	-70.49
+3310.546881	-70.73
+3316.406256	-70.97
+3322.265631	-70.49
+3328.125006	-68.40
+3333.984381	-68.68
+3339.843756	-69.50
+3345.703131	-70.71
+3351.562506	-72.02
+3357.421881	-72.03
+3363.281256	-71.30
+3369.140631	-69.75
+3375.000006	-69.95
+3380.859381	-70.17
+3386.718756	-70.10
+3392.578131	-70.47
+3398.437506	-69.54
+3404.296881	-68.83
+3410.156256	-69.57
+3416.015631	-69.21
+3421.875006	-68.25
+3427.734381	-69.57
+3433.593756	-70.04
+3439.453131	-70.76
+3445.312506	-69.62
+3451.171881	-67.04
+3457.031256	-66.41
+3462.890631	-67.48
+3468.750006	-68.68
+3474.609381	-68.66
+3480.468756	-68.71
+3486.328131	-69.65
+3492.187506	-68.06
+3498.046881	-67.10
+3503.906256	-66.22
+3509.765631	-67.67
+3515.625006	-67.98
+3521.484381	-68.83
+3527.343756	-68.33
+3533.203131	-68.24
+3539.062506	-66.75
+3544.921881	-66.33
+3550.781256	-67.34
+3556.640631	-67.95
+3562.500006	-67.40
+3568.359381	-66.79
+3574.218756	-66.62
+3580.078131	-66.64
+3585.937506	-66.35
+3591.796881	-67.08
+3597.656256	-67.72
+3603.515631	-68.62
+3609.375006	-69.49
+3615.234381	-69.77
+3621.093756	-69.84
+3626.953131	-70.56
+3632.812506	-71.96
+3638.671881	-72.72
+3644.531256	-72.34
+3650.390631	-72.47
+3656.250006	-72.88
+3662.109381	-73.24
+3667.968756	-72.98
+3673.828131	-72.02
+3679.687506	-72.83
+3685.546881	-74.49
+3691.406256	-74.30
+3697.265631	-74.13
+3703.125006	-74.95
+3708.984381	-75.04
+3714.843756	-73.99
+3720.703131	-73.41
+3726.562506	-74.05
+3732.421881	-74.39
+3738.281256	-74.69
+3744.140631	-75.04
+3750.000006	-74.25
+3755.859381	-73.31
+3761.718756	-73.34
+3767.578131	-74.38
+3773.437506	-74.34
+3779.296881	-74.09
+3785.156256	-73.51
+3791.015631	-73.29
+3796.875006	-72.94
+3802.734381	-73.49
+3808.593756	-73.15
+3814.453132	-73.71
+3820.312507	-74.26
+3826.171882	-73.54
+3832.031257	-72.65
+3837.890632	-72.59
+3843.750007	-73.88
+3849.609382	-74.15
+3855.468757	-74.54
+3861.328132	-73.89
+3867.187507	-73.01
+3873.046882	-73.30
+3878.906257	-72.98
+3884.765632	-72.62
+3890.625007	-73.71
+3896.484382	-73.65
+3902.343757	-72.65
+3908.203132	-73.30
+3914.062507	-74.00
+3919.921882	-73.81
+3925.781257	-74.55
+3931.640632	-76.11
+3937.500007	-75.93
+3943.359382	-75.18
+3949.218757	-74.42
+3955.078132	-74.36
+3960.937507	-74.85
+3966.796882	-74.81
+3972.656257	-74.39
+3978.515632	-74.45
+3984.375007	-73.51
+3990.234382	-73.27
+3996.093757	-73.13
+4001.953132	-71.53
+4007.812507	-71.07
+4013.671882	-72.12
+4019.531257	-72.44
+4025.390632	-72.48
+4031.250007	-73.31
+4037.109382	-73.93
+4042.968757	-74.05
+4048.828132	-74.50
+4054.687507	-74.35
+4060.546882	-75.03
+4066.406257	-74.50
+4072.265632	-73.11
+4078.125007	-74.18
+4083.984382	-74.69
+4089.843757	-73.67
+4095.703132	-73.95
+4101.562507	-76.15
+4107.421882	-75.95
+4113.281257	-72.70
+4119.140632	-72.51
+4125.000007	-73.60
+4130.859382	-73.38
+4136.718757	-74.63
+4142.578132	-78.12
+4148.437507	-78.03
+4154.296882	-77.51
+4160.156257	-76.33
+4166.015632	-72.76
+4171.875007	-69.69
+4177.734382	-72.15
+4183.593757	-75.28
+4189.453132	-75.86
+4195.312507	-75.30
+4201.171882	-73.15
+4207.031257	-73.58
+4212.890632	-72.69
+4218.750007	-74.01
+4224.609382	-76.97
+4230.468757	-76.84
+4236.328132	-77.25
+4242.187507	-75.35
+4248.046882	-74.87
+4253.906257	-75.32
+4259.765632	-74.98
+4265.625007	-75.67
+4271.484382	-74.25
+4277.343757	-73.12
+4283.203132	-74.21
+4289.062507	-74.74
+4294.921882	-75.16
+4300.781257	-75.02
+4306.640632	-73.10
+4312.500007	-73.21
+4318.359382	-71.81
+4324.218757	-71.89
+4330.078132	-72.68
+4335.937507	-72.29
+4341.796882	-72.95
+4347.656257	-73.55
+4353.515632	-74.20
+4359.375007	-73.73
+4365.234382	-73.51
+4371.093757	-73.63
+4376.953132	-73.37
+4382.812507	-73.28
+4388.671882	-74.03
+4394.531257	-73.43
+4400.390633	-71.49
+4406.250008	-70.25
+4412.109383	-69.93
+4417.968758	-70.89
+4423.828133	-71.86
+4429.687508	-71.13
+4435.546883	-70.74
+4441.406258	-71.24
+4447.265633	-73.21
+4453.125008	-73.42
+4458.984383	-73.58
+4464.843758	-73.87
+4470.703133	-73.73
+4476.562508	-73.65
+4482.421883	-75.47
+4488.281258	-75.99
+4494.140633	-73.89
+4500.000008	-73.87
+4505.859383	-74.10
+4511.718758	-75.39
+4517.578133	-75.96
+4523.437508	-76.27
+4529.296883	-77.49
+4535.156258	-77.07
+4541.015633	-76.04
+4546.875008	-75.20
+4552.734383	-76.91
+4558.593758	-77.94
+4564.453133	-77.35
+4570.312508	-78.23
+4576.171883	-76.21
+4582.031258	-75.23
+4587.890633	-75.12
+4593.750008	-74.86
+4599.609383	-73.26
+4605.468758	-73.64
+4611.328133	-73.59
+4617.187508	-74.86
+4623.046883	-73.59
+4628.906258	-73.78
+4634.765633	-74.05
+4640.625008	-72.53
+4646.484383	-73.34
+4652.343758	-72.82
+4658.203133	-73.05
+4664.062508	-72.45
+4669.921883	-73.17
+4675.781258	-73.32
+4681.640633	-73.64
+4687.500008	-73.03
+4693.359383	-72.41
+4699.218758	-72.76
+4705.078133	-72.76
+4710.937508	-74.48
+4716.796883	-73.93
+4722.656258	-73.86
+4728.515633	-72.35
+4734.375008	-71.98
+4740.234383	-72.95
+4746.093758	-72.71
+4751.953133	-72.62
+4757.812508	-73.25
+4763.671883	-74.41
+4769.531258	-75.25
+4775.390633	-75.04
+4781.250008	-74.87
+4787.109383	-73.77
+4792.968758	-73.42
+4798.828133	-72.60
+4804.687508	-72.85
+4810.546883	-73.96
+4816.406258	-74.06
+4822.265633	-75.09
+4828.125008	-74.71
+4833.984383	-73.44
+4839.843758	-72.33
+4845.703133	-72.54
+4851.562508	-73.68
+4857.421883	-75.09
+4863.281258	-75.51
+4869.140633	-75.29
+4875.000008	-73.60
+4880.859383	-72.05
+4886.718758	-71.60
+4892.578133	-70.58
+4898.437508	-70.02
+4904.296883	-71.15
+4910.156258	-71.77
+4916.015633	-71.47
+4921.875008	-71.36
+4927.734383	-72.64
+4933.593758	-72.99
+4939.453133	-73.31
+4945.312508	-73.24
+4951.171883	-72.70
+4957.031258	-72.54
+4962.890633	-72.85
+4968.750008	-72.79
+4974.609383	-72.56
+4980.468759	-72.77
+4986.328134	-73.18
+4992.187509	-73.92
+4998.046884	-73.12
+5003.906259	-72.93
+5009.765634	-73.21
+5015.625009	-74.11
+5021.484384	-74.29
+5027.343759	-71.48
+5033.203134	-70.31
+5039.062509	-71.03
+5044.921884	-71.83
+5050.781259	-73.72
+5056.640634	-74.94
+5062.500009	-74.03
+5068.359384	-75.03
+5074.218759	-73.09
+5080.078134	-70.33
+5085.937509	-68.87
+5091.796884	-68.44
+5097.656259	-68.60
+5103.515634	-68.89
+5109.375009	-69.91
+5115.234384	-70.30
+5121.093759	-70.73
+5126.953134	-70.08
+5132.812509	-69.09
+5138.671884	-69.65
+5144.531259	-70.80
+5150.390634	-70.66
+5156.250009	-71.35
+5162.109384	-72.33
+5167.968759	-72.08
+5173.828134	-72.84
+5179.687509	-73.16
+5185.546884	-72.04
+5191.406259	-69.74
+5197.265634	-69.12
+5203.125009	-70.11
+5208.984384	-70.73
+5214.843759	-68.82
+5220.703134	-67.68
+5226.562509	-68.22
+5232.421884	-69.32
+5238.281259	-69.40
+5244.140634	-69.39
+5250.000009	-70.53
+5255.859384	-72.92
+5261.718759	-73.36
+5267.578134	-74.31
+5273.437509	-75.11
+5279.296884	-75.05
+5285.156259	-73.82
+5291.015634	-74.56
+5296.875009	-74.89
+5302.734384	-74.39
+5308.593759	-72.30
+5314.453134	-71.13
+5320.312509	-68.96
+5326.171884	-67.99
+5332.031259	-67.00
+5337.890634	-66.08
+5343.750009	-66.11
+5349.609384	-66.87
+5355.468759	-66.80
+5361.328134	-65.23
+5367.187509	-64.98
+5373.046884	-66.71
+5378.906259	-68.72
+5384.765634	-66.93
+5390.625009	-64.01
+5396.484384	-63.44
+5402.343759	-65.00
+5408.203134	-64.68
+5414.062509	-64.43
+5419.921884	-65.52
+5425.781259	-65.42
+5431.640634	-64.59
+5437.500009	-64.39
+5443.359384	-64.39
+5449.218759	-64.97
+5455.078134	-67.21
+5460.937509	-68.14
+5466.796884	-66.09
+5472.656259	-65.19
+5478.515634	-65.70
+5484.375009	-67.61
+5490.234384	-67.29
+5496.093759	-63.94
+5501.953134	-62.37
+5507.812509	-64.26
+5513.671884	-65.95
+5519.531259	-63.55
+5525.390634	-60.52
+5531.250009	-59.70
+5537.109384	-61.23
+5542.968759	-64.38
+5548.828134	-65.01
+5554.687509	-64.76
+5560.546884	-64.46
+5566.406259	-64.59
+5572.265635	-63.38
+5578.12501	-62.48
+5583.984385	-63.15
+5589.84376	-60.26
+5595.703135	-58.11
+5601.56251	-59.01
+5607.421885	-62.22
+5613.28126	-63.85
+5619.140635	-63.19
+5625.00001	-64.04
+5630.859385	-66.99
+5636.71876	-65.45
+5642.578135	-62.05
+5648.43751	-62.44
+5654.296885	-64.25
+5660.15626	-63.00
+5666.015635	-62.22
+5671.87501	-63.65
+5677.734385	-66.19
+5683.59376	-67.35
+5689.453135	-66.10
+5695.31251	-65.55
+5701.171885	-65.53
+5707.03126	-66.36
+5712.890635	-65.76
+5718.75001	-65.60
+5724.609385	-65.57
+5730.46876	-63.07
+5736.328135	-61.89
+5742.18751	-62.35
+5748.046885	-60.98
+5753.90626	-59.92
+5759.765635	-60.62
+5765.62501	-63.63
+5771.484385	-68.33
+5777.34376	-70.46
+5783.203135	-67.38
+5789.06251	-64.94
+5794.921885	-62.81
+5800.78126	-62.82
+5806.640635	-63.12
+5812.50001	-61.72
+5818.359385	-62.12
+5824.21876	-63.34
+5830.078135	-62.84
+5835.93751	-63.24
+5841.796885	-66.95
+5847.65626	-66.74
+5853.515635	-67.05
+5859.37501	-70.41
+5865.234385	-67.88
+5871.09376	-65.01
+5876.953135	-65.73
+5882.81251	-68.56
+5888.671885	-66.27
+5894.53126	-64.96
+5900.390635	-67.81
+5906.25001	-72.13
+5912.109385	-70.60
+5917.96876	-68.24
+5923.828135	-67.28
+5929.68751	-68.16
+5935.546885	-69.84
+5941.40626	-72.28
+5947.265635	-73.38
+5953.12501	-71.52
+5958.984385	-72.08
+5964.84376	-70.61
+5970.703135	-68.25
+5976.56251	-66.82
+5982.421885	-67.18
+5988.28126	-68.56
+5994.140635	-68.20
+6000.00001	-68.16
+6005.859385	-69.26
+6011.71876	-68.20
+6017.578135	-67.76
+6023.43751	-69.07
+6029.296885	-70.26
+6035.15626	-68.29
+6041.015635	-68.11
+6046.87501	-71.40
+6052.734385	-69.24
+6058.59376	-65.86
+6064.453135	-66.53
+6070.31251	-67.86
+6076.171885	-66.00
+6082.03126	-64.70
+6087.890635	-65.60
+6093.75001	-67.32
+6099.609385	-66.41
+6105.46876	-66.19
+6111.328135	-67.58
+6117.18751	-65.36
+6123.046885	-62.87
+6128.90626	-63.28
+6134.765635	-64.70
+6140.62501	-66.41
+6146.484385	-67.77
+6152.34376	-68.76
+6158.203136	-67.52
+6164.062511	-66.92
+6169.921886	-64.96
+6175.781261	-63.97
+6181.640636	-64.60
+6187.500011	-66.42
+6193.359386	-68.07
+6199.218761	-69.65
+6205.078136	-70.26
+6210.937511	-69.21
+6216.796886	-68.22
+6222.656261	-70.27
+6228.515636	-73.61
+6234.375011	-75.26
+6240.234386	-72.53
+6246.093761	-71.26
+6251.953136	-73.72
+6257.812511	-75.68
+6263.671886	-71.35
+6269.531261	-70.67
+6275.390636	-71.30
+6281.250011	-72.91
+6287.109386	-74.55
+6292.968761	-70.66
+6298.828136	-69.92
+6304.687511	-73.21
+6310.546886	-71.11
+6316.406261	-66.22
+6322.265636	-63.68
+6328.125011	-63.41
+6333.984386	-65.36
+6339.843761	-68.04
+6345.703136	-69.84
+6351.562511	-71.01
+6357.421886	-70.51
+6363.281261	-73.30
+6369.140636	-74.39
+6375.000011	-69.68
+6380.859386	-66.25
+6386.718761	-65.77
+6392.578136	-67.87
+6398.437511	-70.49
+6404.296886	-70.74
+6410.156261	-69.83
+6416.015636	-68.33
+6421.875011	-68.65
+6427.734386	-68.30
+6433.593761	-64.51
+6439.453136	-62.92
+6445.312511	-64.78
+6451.171886	-68.97
+6457.031261	-73.33
+6462.890636	-70.93
+6468.750011	-69.18
+6474.609386	-70.33
+6480.468761	-72.85
+6486.328136	-74.59
+6492.187511	-75.03
+6498.046886	-74.15
+6503.906261	-70.77
+6509.765636	-69.60
+6515.625011	-70.37
+6521.484386	-70.23
+6527.343761	-70.76
+6533.203136	-71.54
+6539.062511	-74.15
+6544.921886	-78.95
+6550.781261	-77.14
+6556.640636	-73.83
+6562.500011	-71.72
+6568.359386	-70.93
+6574.218761	-72.11
+6580.078136	-72.33
+6585.937511	-71.64
+6591.796886	-74.82
+6597.656261	-78.77
+6603.515636	-75.96
+6609.375011	-74.06
+6615.234386	-71.37
+6621.093761	-70.67
+6626.953136	-71.94
+6632.812511	-72.91
+6638.671886	-73.03
+6644.531261	-71.72
+6650.390636	-71.26
+6656.250011	-70.95
+6662.109386	-70.37
+6667.968761	-73.51
+6673.828136	-75.39
+6679.687511	-74.89
+6685.546886	-75.77
+6691.406261	-75.43
+6697.265636	-75.49
+6703.125011	-73.06
+6708.984386	-71.61
+6714.843761	-72.36
+6720.703136	-74.01
+6726.562511	-75.80
+6732.421886	-78.68
+6738.281262	-76.56
+6744.140637	-73.35
+6750.000012	-72.03
+6755.859387	-72.50
+6761.718762	-73.79
+6767.578137	-75.55
+6773.437512	-74.95
+6779.296887	-74.19
+6785.156262	-75.91
+6791.015637	-75.11
+6796.875012	-73.30
+6802.734387	-72.86
+6808.593762	-71.39
+6814.453137	-71.71
+6820.312512	-74.39
+6826.171887	-76.18
+6832.031262	-74.70
+6837.890637	-73.33
+6843.750012	-71.26
+6849.609387	-70.26
+6855.468762	-69.84
+6861.328137	-68.78
+6867.187512	-68.37
+6873.046887	-70.19
+6878.906262	-73.10
+6884.765637	-74.30
+6890.625012	-76.18
+6896.484387	-78.56
+6902.343762	-76.16
+6908.203137	-74.10
+6914.062512	-75.12
+6919.921887	-75.06
+6925.781262	-73.28
+6931.640637	-72.78
+6937.500012	-74.95
+6943.359387	-77.44
+6949.218762	-77.44
+6955.078137	-72.50
+6960.937512	-69.70
+6966.796887	-70.00
+6972.656262	-72.55
+6978.515637	-77.06
+6984.375012	-79.95
+6990.234387	-78.50
+6996.093762	-77.19
+7001.953137	-76.94
+7007.812512	-77.50
+7013.671887	-78.20
+7019.531262	-75.80
+7025.390637	-73.85
+7031.250012	-74.04
+7037.109387	-75.27
+7042.968762	-75.91
+7048.828137	-76.65
+7054.687512	-76.06
+7060.546887	-74.36
+7066.406262	-74.00
+7072.265637	-73.63
+7078.125012	-73.45
+7083.984387	-73.37
+7089.843762	-74.13
+7095.703137	-75.27
+7101.562512	-75.59
+7107.421887	-74.86
+7113.281262	-74.25
+7119.140637	-74.93
+7125.000012	-78.79
+7130.859387	-81.06
+7136.718762	-78.24
+7142.578137	-76.26
+7148.437512	-75.98
+7154.296887	-77.05
+7160.156262	-75.50
+7166.015637	-72.87
+7171.875012	-73.00
+7177.734387	-74.10
+7183.593762	-75.47
+7189.453137	-75.16
+7195.312512	-74.16
+7201.171887	-73.68
+7207.031262	-73.68
+7212.890637	-73.24
+7218.750012	-73.54
+7224.609387	-76.11
+7230.468762	-78.43
+7236.328137	-77.56
+7242.187512	-75.25
+7248.046887	-74.01
+7253.906262	-72.91
+7259.765637	-72.77
+7265.625012	-73.15
+7271.484387	-74.41
+7277.343762	-76.71
+7283.203137	-79.19
+7289.062512	-80.05
+7294.921887	-78.48
+7300.781262	-74.49
+7306.640637	-72.50
+7312.500012	-72.94
+7318.359387	-75.05
+7324.218762	-74.62
+7330.078138	-72.82
+7335.937513	-71.79
+7341.796888	-72.31
+7347.656263	-71.96
+7353.515638	-71.78
+7359.375013	-73.44
+7365.234388	-76.28
+7371.093763	-79.22
+7376.953138	-76.73
+7382.812513	-73.21
+7388.671888	-73.14
+7394.531263	-75.07
+7400.390638	-76.42
+7406.250013	-76.86
+7412.109388	-76.45
+7417.968763	-77.29
+7423.828138	-77.74
+7429.687513	-76.13
+7435.546888	-74.70
+7441.406263	-72.81
+7447.265638	-72.06
+7453.125013	-73.38
+7458.984388	-75.39
+7464.843763	-74.90
+7470.703138	-74.70
+7476.562513	-75.89
+7482.421888	-75.23
+7488.281263	-74.06
+7494.140638	-73.33
+7500.000013	-74.17
+7505.859388	-73.96
+7511.718763	-74.06
+7517.578138	-76.74
+7523.437513	-76.51
+7529.296888	-74.26
+7535.156263	-73.25
+7541.015638	-72.49
+7546.875013	-71.68
+7552.734388	-72.16
+7558.593763	-72.86
+7564.453138	-71.48
+7570.312513	-70.12
+7576.171888	-71.17
+7582.031263	-74.89
+7587.890638	-70.99
+7593.750013	-68.09
+7599.609388	-69.69
+7605.468763	-72.80
+7611.328138	-72.37
+7617.187513	-72.83
+7623.046888	-74.36
+7628.906263	-74.22
+7634.765638	-73.44
+7640.625013	-73.39
+7646.484388	-72.57
+7652.343763	-70.86
+7658.203138	-70.43
+7664.062513	-71.30
+7669.921888	-70.85
+7675.781263	-68.78
+7681.640638	-69.31
+7687.500013	-72.82
+7693.359388	-73.74
+7699.218763	-73.56
+7705.078138	-76.74
+7710.937513	-78.33
+7716.796888	-75.97
+7722.656263	-74.17
+7728.515638	-73.20
+7734.375013	-74.28
+7740.234388	-74.92
+7746.093763	-70.54
+7751.953138	-68.84
+7757.812513	-69.54
+7763.671888	-70.92
+7769.531263	-72.40
+7775.390638	-73.71
+7781.250013	-74.29
+7787.109388	-73.84
+7792.968763	-72.61
+7798.828138	-71.16
+7804.687513	-70.70
+7810.546888	-72.07
+7816.406263	-76.84
+7822.265638	-78.76
+7828.125013	-76.92
+7833.984388	-76.63
+7839.843763	-74.78
+7845.703138	-72.80
+7851.562513	-72.64
+7857.421888	-74.03
+7863.281263	-75.76
+7869.140638	-76.77
+7875.000013	-77.47
+7880.859388	-77.37
+7886.718763	-78.22
+7892.578138	-76.94
+7898.437513	-73.32
+7904.296888	-72.39
+7910.156264	-74.62
+7916.015639	-78.21
+7921.875014	-75.41
+7927.734389	-73.24
+7933.593764	-74.88
+7939.453139	-77.93
+7945.312514	-75.35
+7951.171889	-72.40
+7957.031264	-72.30
+7962.890639	-75.89
+7968.750014	-79.24
+7974.609389	-77.24
+7980.468764	-76.57
+7986.328139	-75.20
+7992.187514	-75.32
+7998.046889	-76.76
+8003.906264	-75.82
+8009.765639	-74.16
+8015.625014	-74.18
+8021.484389	-75.60
+8027.343764	-77.57
+8033.203139	-73.84
+8039.062514	-71.09
+8044.921889	-70.53
+8050.781264	-70.60
+8056.640639	-69.82
+8062.500014	-69.55
+8068.359389	-72.00
+8074.218764	-76.77
+8080.078139	-74.53
+8085.937514	-71.94
+8091.796889	-70.69
+8097.656264	-69.68
+8103.515639	-68.45
+8109.375014	-68.25
+8115.234389	-69.00
+8121.093764	-71.39
+8126.953139	-75.74
+8132.812514	-76.91
+8138.671889	-75.68
+8144.531264	-75.21
+8150.390639	-73.06
+8156.250014	-71.51
+8162.109389	-72.05
+8167.968764	-74.13
+8173.828139	-78.65
+8179.687514	-79.72
+8185.546889	-76.46
+8191.406264	-73.46
+8197.265639	-73.90
+8203.125014	-77.46
+8208.984389	-75.92
+8214.843764	-73.70
+8220.703139	-73.02
+8226.562514	-72.74
+8232.421889	-70.69
+8238.281264	-69.89
+8244.140639	-72.02
+8250.000014	-75.83
+8255.859389	-78.30
+8261.718764	-77.86
+8267.578139	-76.25
+8273.437514	-75.81
+8279.296889	-76.41
+8285.156264	-76.33
+8291.015639	-75.63
+8296.875014	-74.19
+8302.734389	-73.35
+8308.593764	-73.71
+8314.453139	-76.61
+8320.312514	-79.68
+8326.171889	-80.02
+8332.031264	-79.03
+8337.890639	-78.69
+8343.750014	-79.88
+8349.609389	-78.21
+8355.468764	-77.04
+8361.328139	-77.08
+8367.187514	-77.40
+8373.046889	-75.57
+8378.906264	-74.71
+8384.765639	-77.55
+8390.625014	-75.88
+8396.484389	-73.23
+8402.343764	-75.17
+8408.203139	-76.81
+8414.062514	-73.81
+8419.921889	-73.64
+8425.781264	-75.21
+8431.640639	-74.89
+8437.500014	-75.25
+8443.359389	-76.60
+8449.218764	-76.50
+8455.078139	-76.43
+8460.937514	-74.69
+8466.796889	-73.89
+8472.656264	-73.97
+8478.515639	-73.32
+8484.375014	-72.27
+8490.234389	-73.52
+8496.093765	-76.87
+8501.95314	-76.57
+8507.812515	-74.64
+8513.67189	-75.81
+8519.531265	-76.33
+8525.39064	-74.29
+8531.250015	-73.21
+8537.10939	-73.05
+8542.968765	-74.11
+8548.82814	-75.68
+8554.687515	-78.68
+8560.54689	-79.52
+8566.406265	-77.15
+8572.26564	-75.90
+8578.125015	-76.47
+8583.98439	-78.54
+8589.843765	-76.02
+8595.70314	-74.03
+8601.562515	-74.01
+8607.42189	-76.15
+8613.281265	-78.84
+8619.14064	-79.67
+8625.000015	-79.43
+8630.85939	-78.51
+8636.718765	-77.54
+8642.57814	-76.98
+8648.437515	-76.42
+8654.29689	-75.57
+8660.156265	-76.13
+8666.01564	-78.88
+8671.875015	-80.59
+8677.73439	-80.54
+8683.593765	-77.03
+8689.45314	-75.50
+8695.312515	-76.65
+8701.17189	-77.81
+8707.031265	-76.95
+8712.89064	-77.41
+8718.750015	-77.27
+8724.60939	-78.15
+8730.468765	-78.90
+8736.32814	-78.60
+8742.187515	-79.00
+8748.04689	-81.10
+8753.906265	-79.19
+8759.76564	-77.47
+8765.625015	-75.86
+8771.48439	-75.12
+8777.343765	-75.71
+8783.20314	-77.64
+8789.062515	-79.06
+8794.92189	-77.61
+8800.781265	-76.50
+8806.64064	-77.48
+8812.500015	-78.76
+8818.35939	-77.14
+8824.218765	-77.15
+8830.07814	-78.92
+8835.937515	-78.38
+8841.79689	-77.75
+8847.656265	-77.85
+8853.51564	-78.41
+8859.375015	-76.93
+8865.23439	-74.49
+8871.093765	-74.60
+8876.95314	-77.81
+8882.812515	-79.03
+8888.67189	-79.29
+8894.531265	-80.05
+8900.39064	-80.53
+8906.250015	-80.44
+8912.10939	-78.40
+8917.968765	-77.14
+8923.82814	-78.01
+8929.687515	-78.19
+8935.54689	-79.77
+8941.406265	-81.49
+8947.26564	-79.61
+8953.125015	-77.01
+8958.98439	-75.98
+8964.843765	-77.06
+8970.70314	-78.46
+8976.562515	-78.48
+8982.42189	-79.52
+8988.281265	-79.20
+8994.14064	-77.82
+9000.000015	-78.06
+9005.85939	-77.00
+9011.718765	-75.44
+9017.57814	-76.84
+9023.437515	-78.75
+9029.29689	-79.69
+9035.156265	-80.83
+9041.01564	-79.14
+9046.875015	-78.56
+9052.73439	-77.88
+9058.593765	-78.37
+9064.45314	-78.29
+9070.312515	-77.23
+9076.17189	-77.47
+9082.031266	-77.60
+9087.890641	-77.86
+9093.750016	-79.28
+9099.609391	-78.51
+9105.468766	-77.52
+9111.328141	-78.29
+9117.187516	-78.20
+9123.046891	-77.95
+9128.906266	-77.87
+9134.765641	-78.50
+9140.625016	-79.98
+9146.484391	-79.40
+9152.343766	-76.80
+9158.203141	-77.06
+9164.062516	-77.15
+9169.921891	-76.51
+9175.781266	-76.48
+9181.640641	-77.08
+9187.500016	-78.62
+9193.359391	-78.12
+9199.218766	-76.74
+9205.078141	-76.55
+9210.937516	-77.90
+9216.796891	-77.78
+9222.656266	-76.11
+9228.515641	-75.29
+9234.375016	-73.24
+9240.234391	-73.31
+9246.093766	-75.88
+9251.953141	-77.50
+9257.812516	-78.33
+9263.671891	-78.40
+9269.531266	-77.08
+9275.390641	-74.86
+9281.250016	-73.76
+9287.109391	-75.71
+9292.968766	-77.65
+9298.828141	-76.72
+9304.687516	-76.79
+9310.546891	-78.08
+9316.406266	-78.13
+9322.265641	-78.75
+9328.125016	-78.75
+9333.984391	-75.26
+9339.843766	-73.74
+9345.703141	-75.38
+9351.562516	-76.01
+9357.421891	-74.25
+9363.281266	-74.74
+9369.140641	-76.46
+9375.000016	-77.03
+9380.859391	-76.02
+9386.718766	-76.19
+9392.578141	-78.24
+9398.437516	-76.81
+9404.296891	-76.56
+9410.156266	-76.19
+9416.015641	-75.39
+9421.875016	-76.76
+9427.734391	-78.70
+9433.593766	-78.39
+9439.453141	-76.86
+9445.312516	-76.38
+9451.171891	-76.19
+9457.031266	-77.13
+9462.890641	-76.37
+9468.750016	-74.04
+9474.609391	-71.24
+9480.468766	-70.39
+9486.328141	-72.79
+9492.187516	-75.73
+9498.046891	-74.43
+9503.906266	-75.19
+9509.765641	-77.27
+9515.625016	-77.47
+9521.484391	-77.75
+9527.343766	-77.99
+9533.203141	-73.50
+9539.062516	-69.97
+9544.921891	-69.91
+9550.781266	-73.10
+9556.640641	-77.72
+9562.500016	-79.52
+9568.359391	-79.20
+9574.218766	-77.13
+9580.078141	-75.59
+9585.937516	-71.93
+9591.796891	-69.79
+9597.656266	-70.63
+9603.515641	-74.69
+9609.375016	-80.15
+9615.234391	-81.21
+9621.093766	-80.19
+9626.953141	-78.43
+9632.812516	-76.90
+9638.671891	-74.16
+9644.531266	-73.10
+9650.390641	-74.90
+9656.250016	-77.34
+9662.109391	-76.84
+9667.968767	-77.59
+9673.828142	-79.46
+9679.687517	-77.02
+9685.546892	-76.45
+9691.406267	-79.24
+9697.265642	-81.42
+9703.125017	-79.12
+9708.984392	-78.03
+9714.843767	-76.96
+9720.703142	-76.42
+9726.562517	-73.03
+9732.421892	-71.52
+9738.281267	-73.41
+9744.140642	-77.42
+9750.000017	-74.87
+9755.859392	-73.32
+9761.718767	-74.74
+9767.578142	-78.33
+9773.437517	-79.38
+9779.296892	-76.11
+9785.156267	-74.21
+9791.015642	-72.55
+9796.875017	-72.13
+9802.734392	-72.08
+9808.593767	-72.86
+9814.453142	-73.58
+9820.312517	-74.12
+9826.171892	-76.17
+9832.031267	-79.20
+9837.890642	-77.78
+9843.750017	-77.34
+9849.609392	-79.06
+9855.468767	-77.43
+9861.328142	-75.16
+9867.187517	-74.16
+9873.046892	-75.07
+9878.906267	-77.14
+9884.765642	-77.49
+9890.625017	-74.73
+9896.484392	-72.39
+9902.343767	-72.05
+9908.203142	-73.16
+9914.062517	-75.99
+9919.921892	-76.72
+9925.781267	-75.98
+9931.640642	-76.34
+9937.500017	-78.20
+9943.359392	-79.38
+9949.218767	-77.51
+9955.078142	-75.58
+9960.937517	-76.62
+9966.796892	-77.68
+9972.656267	-75.17
+9978.515642	-74.65
+9984.375017	-76.74
+9990.234392	-76.89
+9996.093767	-76.47
+10001.953142	-75.07
+10007.812517	-75.26
+10013.671892	-76.67
+10019.531267	-78.11
+10025.390642	-76.58
+10031.250017	-76.04
+10037.109392	-77.76
+10042.968767	-79.48
+10048.828142	-80.12
+10054.687517	-82.27
+10060.546892	-82.45
+10066.406267	-80.57
+10072.265642	-81.44
+10078.125017	-80.11
+10083.984392	-75.94
+10089.843767	-74.12
+10095.703142	-74.84
+10101.562517	-77.63
+10107.421892	-78.78
+10113.281267	-79.43
+10119.140642	-79.84
+10125.000017	-79.76
+10130.859392	-79.32
+10136.718767	-78.30
+10142.578142	-77.50
+10148.437517	-78.10
+10154.296892	-79.10
+10160.156267	-80.64
+10166.015642	-82.81
+10171.875017	-82.58
+10177.734392	-80.47
+10183.593767	-76.50
+10189.453142	-73.44
+10195.312517	-73.36
+10201.171892	-75.67
+10207.031267	-77.68
+10212.890642	-78.25
+10218.750017	-79.99
+10224.609392	-80.89
+10230.468767	-78.17
+10236.328142	-76.58
+10242.187517	-78.05
+10248.046892	-78.67
+10253.906268	-75.72
+10259.765643	-75.50
+10265.625018	-77.81
+10271.484393	-76.64
+10277.343768	-75.53
+10283.203143	-77.32
+10289.062518	-79.47
+10294.921893	-76.95
+10300.781268	-76.40
+10306.640643	-79.29
+10312.500018	-77.42
+10318.359393	-73.37
+10324.218768	-72.62
+10330.078143	-74.30
+10335.937518	-77.62
+10341.796893	-80.48
+10347.656268	-80.67
+10353.515643	-78.48
+10359.375018	-75.86
+10365.234393	-74.96
+10371.093768	-76.61
+10376.953143	-79.74
+10382.812518	-82.36
+10388.671893	-79.20
+10394.531268	-77.27
+10400.390643	-78.39
+10406.250018	-78.18
+10412.109393	-77.90
+10417.968768	-77.26
+10423.828143	-76.88
+10429.687518	-76.78
+10435.546893	-76.32
+10441.406268	-76.25
+10447.265643	-76.20
+10453.125018	-77.65
+10458.984393	-79.16
+10464.843768	-79.37
+10470.703143	-78.03
+10476.562518	-78.52
+10482.421893	-79.35
+10488.281268	-78.97
+10494.140643	-79.43
+10500.000018	-77.44
+10505.859393	-77.37
+10511.718768	-77.77
+10517.578143	-79.65
+10523.437518	-82.38
+10529.296893	-80.76
+10535.156268	-79.52
+10541.015643	-78.55
+10546.875018	-78.36
+10552.734393	-78.05
+10558.593768	-78.11
+10564.453143	-78.39
+10570.312518	-77.45
+10576.171893	-77.15
+10582.031268	-80.83
+10587.890643	-82.15
+10593.750018	-79.70
+10599.609393	-80.53
+10605.468768	-78.20
+10611.328143	-76.97
+10617.187518	-77.95
+10623.046893	-79.45
+10628.906268	-81.62
+10634.765643	-78.46
+10640.625018	-76.86
+10646.484393	-78.02
+10652.343768	-78.65
+10658.203143	-80.15
+10664.062518	-80.59
+10669.921893	-81.07
+10675.781268	-79.07
+10681.640643	-79.92
+10687.500018	-81.31
+10693.359393	-79.19
+10699.218768	-79.00
+10705.078143	-79.23
+10710.937518	-78.63
+10716.796893	-79.32
+10722.656268	-80.45
+10728.515643	-80.77
+10734.375018	-81.43
+10740.234393	-81.31
+10746.093768	-80.72
+10751.953143	-80.22
+10757.812518	-80.51
+10763.671893	-80.01
+10769.531268	-78.17
+10775.390643	-79.15
+10781.250018	-82.03
+10787.109393	-81.88
+10792.968768	-81.44
+10798.828143	-81.49
+10804.687518	-82.19
+10810.546893	-83.98
+10816.406268	-83.63
+10822.265643	-80.67
+10828.125018	-78.76
+10833.984393	-78.19
+10839.843769	-79.38
+10845.703144	-79.65
+10851.562519	-77.22
+10857.421894	-76.62
+10863.281269	-77.63
+10869.140644	-80.02
+10875.000019	-79.00
+10880.859394	-77.19
+10886.718769	-77.99
+10892.578144	-80.19
+10898.437519	-80.50
+10904.296894	-79.43
+10910.156269	-78.38
+10916.015644	-78.57
+10921.875019	-77.78
+10927.734394	-75.93
+10933.593769	-75.35
+10939.453144	-76.36
+10945.312519	-78.61
+10951.171894	-79.99
+10957.031269	-79.58
+10962.890644	-78.36
+10968.750019	-77.89
+10974.609394	-78.42
+10980.468769	-79.41
+10986.328144	-81.05
+10992.187519	-82.58
+10998.046894	-81.01
+11003.906269	-79.71
+11009.765644	-80.38
+11015.625019	-82.59
+11021.484394	-83.70
+11027.343769	-83.87
+11033.203144	-82.08
+11039.062519	-80.09
+11044.921894	-80.28
+11050.781269	-79.94
+11056.640644	-79.08
+11062.500019	-78.24
+11068.359394	-78.45
+11074.218769	-79.51
+11080.078144	-81.94
+11085.937519	-81.26
+11091.796894	-78.21
+11097.656269	-77.49
+11103.515644	-78.44
+11109.375019	-78.94
+11115.234394	-78.47
+11121.093769	-80.46
+11126.953144	-81.83
+11132.812519	-80.29
+11138.671894	-78.28
+11144.531269	-78.23
+11150.390644	-80.17
+11156.250019	-80.76
+11162.109394	-81.05
+11167.968769	-80.92
+11173.828144	-79.86
+11179.687519	-78.46
+11185.546894	-78.66
+11191.406269	-78.75
+11197.265644	-80.35
+11203.125019	-80.28
+11208.984394	-80.06
+11214.843769	-80.52
+11220.703144	-79.56
+11226.562519	-79.68
+11232.421894	-81.38
+11238.281269	-82.66
+11244.140644	-81.46
+11250.000019	-80.21
+11255.859394	-79.31
+11261.718769	-79.53
+11267.578144	-79.20
+11273.437519	-79.69
+11279.296894	-81.13
+11285.156269	-80.13
+11291.015644	-80.24
+11296.875019	-81.74
+11302.734394	-82.59
+11308.593769	-80.70
+11314.453144	-78.87
+11320.312519	-79.92
+11326.171894	-84.01
+11332.031269	-82.30
+11337.890644	-78.38
+11343.750019	-75.41
+11349.609394	-74.04
+11355.468769	-74.52
+11361.328144	-77.67
+11367.187519	-81.22
+11373.046894	-81.17
+11378.906269	-80.34
+11384.765644	-78.94
+11390.625019	-80.74
+11396.484394	-81.46
+11402.343769	-81.41
+11408.203144	-81.66
+11414.062519	-81.61
+11419.921894	-82.70
+11425.78127	-82.19
+11431.640645	-82.10
+11437.50002	-83.41
+11443.359395	-83.27
+11449.21877	-84.08
+11455.078145	-86.14
+11460.93752	-85.53
+11466.796895	-84.23
+11472.65627	-82.34
+11478.515645	-82.53
+11484.37502	-81.75
+11490.234395	-82.01
+11496.09377	-83.44
+11501.953145	-82.10
+11507.81252	-82.86
+11513.671895	-83.27
+11519.53127	-81.81
+11525.390645	-82.17
+11531.25002	-82.15
+11537.109395	-80.59
+11542.96877	-79.96
+11548.828145	-81.76
+11554.68752	-81.93
+11560.546895	-81.28
+11566.40627	-80.68
+11572.265645	-77.85
+11578.12502	-77.33
+11583.984395	-79.07
+11589.84377	-79.92
+11595.703145	-79.65
+11601.56252	-82.34
+11607.421895	-84.14
+11613.28127	-82.91
+11619.140645	-83.82
+11625.00002	-80.65
+11630.859395	-78.93
+11636.71877	-79.06
+11642.578145	-78.07
+11648.43752	-79.14
+11654.296895	-80.84
+11660.15627	-81.04
+11666.015645	-79.98
+11671.87502	-79.78
+11677.734395	-79.83
+11683.59377	-76.68
+11689.453145	-76.17
+11695.31252	-78.03
+11701.171895	-78.60
+11707.03127	-77.73
+11712.890645	-78.15
+11718.75002	-80.16
+11724.609395	-81.49
+11730.46877	-83.14
+11736.328145	-84.78
+11742.18752	-80.50
+11748.046895	-78.90
+11753.90627	-81.66
+11759.765645	-81.94
+11765.62502	-79.51
+11771.484395	-80.63
+11777.34377	-83.59
+11783.203145	-85.62
+11789.06252	-82.90
+11794.921895	-80.77
+11800.78127	-80.46
+11806.640645	-81.99
+11812.50002	-84.68
+11818.359395	-83.53
+11824.21877	-80.24
+11830.078145	-77.78
+11835.93752	-77.97
+11841.796895	-81.72
+11847.65627	-84.52
+11853.515645	-82.31
+11859.37502	-81.72
+11865.234395	-82.19
+11871.09377	-82.34
+11876.953145	-83.77
+11882.81252	-82.83
+11888.671895	-81.38
+11894.53127	-82.25
+11900.390645	-82.31
+11906.25002	-84.12
+11912.109395	-85.38
+11917.96877	-82.59
+11923.828145	-81.43
+11929.68752	-82.52
+11935.546895	-82.97
+11941.40627	-81.79
+11947.265645	-81.81
+11953.12502	-80.76
+11958.984395	-78.39
+11964.84377	-78.28
+11970.703145	-80.63
+11976.56252	-82.04
+11982.421895	-81.25
+11988.28127	-82.39
+11994.140645	-84.00
+12000.00002	-83.66
+12005.859395	-85.21
+12011.718771	-86.27
+12017.578146	-84.64
+12023.437521	-86.60
+12029.296896	-87.13
+12035.156271	-84.05
+12041.015646	-84.52
+12046.875021	-84.54
+12052.734396	-80.94
+12058.593771	-79.18
+12064.453146	-79.93
+12070.312521	-80.67
+12076.171896	-81.85
+12082.031271	-82.17
+12087.890646	-79.79
+12093.750021	-78.37
+12099.609396	-79.82
+12105.468771	-82.84
+12111.328146	-85.94
+12117.187521	-84.46
+12123.046896	-82.86
+12128.906271	-83.52
+12134.765646	-86.02
+12140.625021	-85.68
+12146.484396	-85.06
+12152.343771	-83.86
+12158.203146	-82.80
+12164.062521	-81.85
+12169.921896	-83.22
+12175.781271	-84.31
+12181.640646	-85.03
+12187.500021	-85.43
+12193.359396	-84.87
+12199.218771	-85.36
+12205.078146	-85.82
+12210.937521	-85.47
+12216.796896	-84.11
+12222.656271	-84.17
+12228.515646	-85.71
+12234.375021	-85.82
+12240.234396	-85.12
+12246.093771	-84.70
+12251.953146	-84.51
+12257.812521	-83.90
+12263.671896	-85.08
+12269.531271	-84.35
+12275.390646	-84.96
+12281.250021	-85.70
+12287.109396	-84.86
+12292.968771	-83.85
+12298.828146	-83.38
+12304.687521	-84.95
+12310.546896	-86.93
+12316.406271	-87.77
+12322.265646	-86.25
+12328.125021	-86.40
+12333.984396	-87.84
+12339.843771	-86.81
+12345.703146	-85.84
+12351.562521	-83.30
+12357.421896	-82.98
+12363.281271	-85.18
+12369.140646	-85.17
+12375.000021	-85.90
+12380.859396	-86.92
+12386.718771	-87.22
+12392.578146	-85.93
+12398.437521	-85.69
+12404.296896	-85.35
+12410.156271	-87.78
+12416.015646	-89.05
+12421.875021	-89.02
+12427.734396	-89.30
+12433.593771	-87.51
+12439.453146	-84.84
+12445.312521	-85.07
+12451.171896	-86.32
+12457.031271	-87.43
+12462.890646	-87.53
+12468.750021	-86.97
+12474.609396	-86.12
+12480.468771	-84.69
+12486.328146	-84.05
+12492.187521	-85.52
+12498.046896	-85.34
+12503.906271	-86.05
+12509.765646	-85.67
+12515.625021	-84.14
+12521.484396	-84.57
+12527.343771	-83.92
+12533.203146	-83.65
+12539.062521	-85.91
+12544.921896	-87.46
+12550.781271	-86.95
+12556.640646	-88.42
+12562.500021	-89.65
+12568.359396	-86.22
+12574.218771	-84.86
+12580.078146	-87.88
+12585.937521	-85.81
+12591.796896	-83.59
+12597.656272	-85.25
+12603.515647	-87.69
+12609.375022	-88.00
+12615.234397	-87.32
+12621.093772	-87.39
+12626.953147	-85.22
+12632.812522	-84.61
+12638.671897	-85.57
+12644.531272	-84.69
+12650.390647	-85.48
+12656.250022	-89.32
+12662.109397	-89.97
+12667.968772	-88.07
+12673.828147	-88.82
+12679.687522	-89.59
+12685.546897	-88.68
+12691.406272	-85.83
+12697.265647	-84.65
+12703.125022	-84.84
+12708.984397	-85.58
+12714.843772	-85.69
+12720.703147	-86.04
+12726.562522	-88.33
+12732.421897	-87.94
+12738.281272	-87.81
+12744.140647	-90.14
+12750.000022	-89.51
+12755.859397	-88.47
+12761.718772	-87.04
+12767.578147	-85.47
+12773.437522	-85.16
+12779.296897	-86.87
+12785.156272	-89.79
+12791.015647	-89.75
+12796.875022	-87.79
+12802.734397	-88.71
+12808.593772	-90.37
+12814.453147	-91.52
+12820.312522	-89.88
+12826.171897	-87.14
+12832.031272	-87.56
+12837.890647	-89.94
+12843.750022	-89.86
+12849.609397	-88.65
+12855.468772	-89.28
+12861.328147	-90.63
+12867.187522	-88.43
+12873.046897	-87.61
+12878.906272	-88.82
+12884.765647	-87.91
+12890.625022	-87.04
+12896.484397	-88.16
+12902.343772	-88.89
+12908.203147	-90.21
+12914.062522	-90.12
+12919.921897	-88.82
+12925.781272	-89.05
+12931.640647	-88.57
+12937.500022	-88.34
+12943.359397	-88.41
+12949.218772	-87.00
+12955.078147	-85.94
+12960.937522	-85.69
+12966.796897	-85.34
+12972.656272	-86.96
+12978.515647	-89.13
+12984.375022	-89.86
+12990.234397	-89.62
+12996.093772	-87.65
+13001.953147	-85.73
+13007.812522	-85.40
+13013.671897	-85.39
+13019.531272	-85.81
+13025.390647	-87.67
+13031.250022	-89.20
+13037.109397	-91.12
+13042.968772	-89.10
+13048.828147	-89.07
+13054.687522	-90.53
+13060.546897	-89.45
+13066.406272	-87.95
+13072.265647	-87.06
+13078.125022	-86.07
+13083.984397	-86.13
+13089.843772	-88.16
+13095.703147	-88.83
+13101.562522	-89.00
+13107.421897	-89.72
+13113.281272	-88.40
+13119.140647	-87.92
+13125.000022	-87.50
+13130.859397	-86.91
+13136.718772	-86.56
+13142.578147	-87.75
+13148.437522	-89.72
+13154.296897	-87.41
+13160.156272	-84.89
+13166.015647	-84.22
+13171.875022	-85.28
+13177.734397	-86.71
+13183.593773	-86.82
+13189.453148	-86.34
+13195.312523	-86.69
+13201.171898	-87.13
+13207.031273	-84.74
+13212.890648	-86.43
+13218.750023	-91.35
+13224.609398	-87.70
+13230.468773	-85.08
+13236.328148	-85.68
+13242.187523	-85.16
+13248.046898	-85.96
+13253.906273	-89.06
+13259.765648	-87.27
+13265.625023	-85.41
+13271.484398	-86.42
+13277.343773	-88.99
+13283.203148	-87.66
+13289.062523	-85.99
+13294.921898	-86.84
+13300.781273	-87.98
+13306.640648	-88.41
+13312.500023	-88.35
+13318.359398	-88.16
+13324.218773	-89.02
+13330.078148	-89.56
+13335.937523	-88.94
+13341.796898	-88.21
+13347.656273	-88.59
+13353.515648	-89.40
+13359.375023	-88.68
+13365.234398	-89.56
+13371.093773	-87.08
+13376.953148	-86.75
+13382.812523	-89.39
+13388.671898	-90.39
+13394.531273	-89.71
+13400.390648	-88.69
+13406.250023	-87.08
+13412.109398	-86.95
+13417.968773	-87.82
+13423.828148	-87.76
+13429.687523	-84.31
+13435.546898	-82.45
+13441.406273	-84.84
+13447.265648	-86.33
+13453.125023	-87.28
+13458.984398	-89.34
+13464.843773	-87.51
+13470.703148	-88.14
+13476.562523	-89.83
+13482.421898	-87.52
+13488.281273	-87.96
+13494.140648	-89.15
+13500.000023	-88.59
+13505.859398	-89.53
+13511.718773	-90.58
+13517.578148	-87.75
+13523.437523	-87.23
+13529.296898	-89.44
+13535.156273	-89.57
+13541.015648	-87.61
+13546.875023	-85.20
+13552.734398	-85.09
+13558.593773	-88.71
+13564.453148	-88.56
+13570.312523	-86.11
+13576.171898	-86.30
+13582.031273	-88.77
+13587.890648	-88.98
+13593.750023	-87.11
+13599.609398	-87.94
+13605.468773	-88.89
+13611.328148	-89.55
+13617.187523	-91.81
+13623.046898	-91.34
+13628.906273	-92.30
+13634.765648	-92.04
+13640.625023	-89.73
+13646.484398	-89.24
+13652.343773	-91.16
+13658.203148	-91.22
+13664.062523	-91.40
+13669.921898	-92.25
+13675.781273	-90.18
+13681.640648	-87.70
+13687.500023	-88.85
+13693.359398	-92.34
+13699.218773	-91.19
+13705.078148	-88.20
+13710.937523	-88.00
+13716.796898	-89.85
+13722.656273	-87.78
+13728.515648	-86.55
+13734.375023	-89.16
+13740.234398	-89.80
+13746.093773	-87.24
+13751.953148	-87.15
+13757.812523	-87.03
+13763.671898	-87.42
+13769.531274	-88.06
+13775.390649	-88.67
+13781.250024	-89.51
+13787.109399	-88.54
+13792.968774	-87.73
+13798.828149	-87.40
+13804.687524	-87.39
+13810.546899	-88.17
+13816.406274	-90.16
+13822.265649	-88.45
+13828.125024	-85.82
+13833.984399	-86.38
+13839.843774	-88.97
+13845.703149	-91.50
+13851.562524	-87.67
+13857.421899	-84.70
+13863.281274	-85.02
+13869.140649	-86.25
+13875.000024	-88.57
+13880.859399	-89.96
+13886.718774	-89.59
+13892.578149	-87.06
+13898.437524	-86.12
+13904.296899	-88.70
+13910.156274	-92.95
+13916.015649	-91.20
+13921.875024	-89.98
+13927.734399	-91.05
+13933.593774	-92.06
+13939.453149	-88.97
+13945.312524	-86.81
+13951.171899	-87.99
+13957.031274	-88.90
+13962.890649	-88.37
+13968.750024	-90.17
+13974.609399	-90.05
+13980.468774	-87.84
+13986.328149	-86.01
+13992.187524	-84.78
+13998.046899	-85.10
+14003.906274	-87.83
+14009.765649	-89.65
+14015.625024	-88.93
+14021.484399	-89.23
+14027.343774	-87.99
+14033.203149	-88.00
+14039.062524	-87.12
+14044.921899	-84.88
+14050.781274	-86.07
+14056.640649	-88.04
+14062.500024	-86.20
+14068.359399	-85.10
+14074.218774	-85.01
+14080.078149	-84.79
+14085.937524	-85.48
+14091.796899	-85.76
+14097.656274	-84.95
+14103.515649	-85.76
+14109.375024	-88.14
+14115.234399	-90.96
+14121.093774	-89.88
+14126.953149	-91.13
+14132.812524	-90.20
+14138.671899	-87.40
+14144.531274	-85.44
+14150.390649	-85.86
+14156.250024	-87.33
+14162.109399	-88.18
+14167.968774	-87.56
+14173.828149	-87.62
+14179.687524	-89.15
+14185.546899	-88.28
+14191.406274	-88.85
+14197.265649	-90.79
+14203.125024	-89.18
+14208.984399	-87.92
+14214.843774	-88.79
+14220.703149	-90.11
+14226.562524	-89.88
+14232.421899	-88.50
+14238.281274	-88.14
+14244.140649	-89.39
+14250.000024	-91.13
+14255.859399	-88.71
+14261.718774	-86.70
+14267.578149	-86.56
+14273.437524	-87.51
+14279.296899	-88.91
+14285.156274	-88.78
+14291.015649	-86.65
+14296.875024	-86.45
+14302.734399	-86.20
+14308.593774	-86.14
+14314.453149	-89.13
+14320.312524	-89.73
+14326.171899	-87.58
+14332.031274	-89.57
+14337.890649	-90.91
+14343.750024	-86.39
+14349.609399	-84.74
+14355.468775	-86.24
+14361.32815	-89.37
+14367.187525	-90.72
+14373.0469	-89.70
+14378.906275	-88.85
+14384.76565	-88.92
+14390.625025	-89.64
+14396.4844	-90.49
+14402.343775	-90.70
+14408.20315	-90.86
+14414.062525	-90.50
+14419.9219	-88.19
+14425.781275	-86.88
+14431.64065	-85.34
+14437.500025	-84.19
+14443.3594	-85.13
+14449.218775	-86.87
+14455.07815	-88.24
+14460.937525	-90.45
+14466.7969	-90.98
+14472.656275	-88.92
+14478.51565	-86.68
+14484.375025	-86.03
+14490.2344	-85.93
+14496.093775	-86.45
+14501.95315	-88.02
+14507.812525	-89.20
+14513.6719	-89.57
+14519.531275	-88.45
+14525.39065	-88.46
+14531.250025	-87.93
+14537.1094	-86.37
+14542.968775	-85.95
+14548.82815	-87.16
+14554.687525	-88.07
+14560.5469	-88.22
+14566.406275	-88.57
+14572.26565	-86.69
+14578.125025	-84.86
+14583.9844	-84.65
+14589.843775	-86.39
+14595.70315	-86.59
+14601.562525	-86.11
+14607.4219	-89.60
+14613.281275	-92.30
+14619.14065	-91.32
+14625.000025	-91.70
+14630.8594	-90.03
+14636.718775	-89.88
+14642.57815	-91.30
+14648.437525	-90.66
+14654.2969	-89.89
+14660.156275	-89.33
+14666.01565	-88.56
+14671.875025	-89.02
+14677.7344	-90.07
+14683.593775	-89.66
+14689.45315	-87.82
+14695.312525	-87.14
+14701.1719	-88.11
+14707.031275	-90.10
+14712.89065	-89.12
+14718.750025	-88.66
+14724.6094	-89.68
+14730.468775	-88.10
+14736.32815	-88.03
+14742.187525	-89.72
+14748.0469	-90.54
+14753.906275	-88.80
+14759.76565	-88.41
+14765.625025	-90.60
+14771.4844	-90.95
+14777.343775	-89.06
+14783.20315	-89.39
+14789.062525	-90.42
+14794.9219	-88.41
+14800.781275	-86.85
+14806.64065	-86.66
+14812.500025	-88.88
+14818.3594	-87.68
+14824.218775	-86.74
+14830.07815	-88.97
+14835.937525	-88.22
+14841.7969	-86.69
+14847.656275	-87.62
+14853.51565	-90.28
+14859.375025	-92.11
+14865.2344	-92.42
+14871.093775	-91.66
+14876.95315	-90.97
+14882.812525	-90.69
+14888.6719	-89.34
+14894.531275	-89.27
+14900.39065	-89.65
+14906.250025	-88.72
+14912.1094	-87.55
+14917.968775	-86.28
+14923.82815	-85.54
+14929.687525	-87.22
+14935.5469	-90.04
+14941.406276	-88.05
+14947.265651	-86.73
+14953.125026	-88.23
+14958.984401	-89.82
+14964.843776	-85.96
+14970.703151	-83.85
+14976.562526	-85.70
+14982.421901	-87.22
+14988.281276	-85.30
+14994.140651	-86.85
+15000.000026	-90.54
+15005.859401	-89.59
+15011.718776	-89.18
+15017.578151	-88.25
+15023.437526	-87.88
+15029.296901	-89.05
+15035.156276	-85.32
+15041.015651	-84.79
+15046.875026	-87.28
+15052.734401	-86.04
+15058.593776	-87.33
+15064.453151	-89.48
+15070.312526	-87.64
+15076.171901	-88.38
+15082.031276	-86.73
+15087.890651	-84.39
+15093.750026	-83.98
+15099.609401	-83.53
+15105.468776	-83.99
+15111.328151	-86.72
+15117.187526	-88.68
+15123.046901	-87.80
+15128.906276	-85.92
+15134.765651	-85.90
+15140.625026	-86.71
+15146.484401	-85.25
+15152.343776	-84.81
+15158.203151	-87.88
+15164.062526	-89.55
+15169.921901	-88.41
+15175.781276	-85.00
+15181.640651	-82.26
+15187.500026	-82.58
+15193.359401	-85.50
+15199.218776	-88.39
+15205.078151	-88.65
+15210.937526	-88.77
+15216.796901	-87.23
+15222.656276	-87.91
+15228.515651	-89.96
+15234.375026	-90.38
+15240.234401	-92.66
+15246.093776	-90.82
+15251.953151	-88.02
+15257.812526	-85.87
+15263.671901	-85.44
+15269.531276	-85.49
+15275.390651	-86.34
+15281.250026	-88.35
+15287.109401	-88.94
+15292.968776	-87.83
+15298.828151	-85.40
+15304.687526	-84.23
+15310.546901	-85.79
+15316.406276	-86.52
+15322.265651	-86.68
+15328.125026	-89.90
+15333.984401	-89.86
+15339.843776	-88.32
+15345.703151	-87.40
+15351.562526	-84.52
+15357.421901	-83.28
+15363.281276	-83.57
+15369.140651	-85.62
+15375.000026	-88.05
+15380.859401	-87.67
+15386.718776	-86.73
+15392.578151	-86.69
+15398.437526	-88.08
+15404.296901	-87.53
+15410.156276	-88.27
+15416.015651	-87.68
+15421.875026	-86.41
+15427.734401	-88.07
+15433.593776	-90.05
+15439.453151	-87.70
+15445.312526	-86.51
+15451.171901	-87.99
+15457.031276	-90.90
+15462.890651	-90.13
+15468.750026	-87.14
+15474.609401	-86.41
+15480.468776	-88.23
+15486.328151	-89.78
+15492.187526	-89.29
+15498.046901	-88.66
+15503.906276	-87.96
+15509.765651	-86.61
+15515.625026	-82.71
+15521.484401	-80.66
+15527.343777	-82.62
+15533.203152	-87.92
+15539.062527	-88.48
+15544.921902	-89.28
+15550.781277	-90.26
+15556.640652	-89.03
+15562.500027	-89.73
+15568.359402	-90.09
+15574.218777	-86.78
+15580.078152	-86.48
+15585.937527	-87.40
+15591.796902	-86.77
+15597.656277	-86.83
+15603.515652	-84.85
+15609.375027	-84.07
+15615.234402	-86.83
+15621.093777	-89.59
+15626.953152	-88.49
+15632.812527	-87.81
+15638.671902	-86.07
+15644.531277	-84.55
+15650.390652	-85.36
+15656.250027	-88.85
+15662.109402	-90.92
+15667.968777	-90.36
+15673.828152	-88.30
+15679.687527	-87.61
+15685.546902	-89.56
+15691.406277	-89.56
+15697.265652	-88.75
+15703.125027	-89.27
+15708.984402	-88.88
+15714.843777	-88.66
+15720.703152	-85.65
+15726.562527	-83.52
+15732.421902	-84.68
+15738.281277	-88.55
+15744.140652	-91.37
+15750.000027	-92.34
+15755.859402	-91.22
+15761.718777	-91.64
+15767.578152	-91.91
+15773.437527	-90.61
+15779.296902	-87.72
+15785.156277	-86.60
+15791.015652	-87.44
+15796.875027	-88.60
+15802.734402	-87.86
+15808.593777	-86.29
+15814.453152	-85.40
+15820.312527	-86.72
+15826.171902	-89.12
+15832.031277	-91.87
+15837.890652	-89.70
+15843.750027	-86.80
+15849.609402	-86.76
+15855.468777	-87.09
+15861.328152	-86.64
+15867.187527	-86.36
+15873.046902	-86.42
+15878.906277	-88.58
+15884.765652	-90.62
+15890.625027	-89.64
+15896.484402	-89.13
+15902.343777	-88.71
+15908.203152	-88.52
+15914.062527	-89.17
+15919.921902	-89.03
+15925.781277	-87.85
+15931.640652	-88.92
+15937.500027	-90.59
+15943.359402	-88.02
+15949.218777	-86.59
+15955.078152	-88.48
+15960.937527	-91.12
+15966.796902	-89.74
+15972.656277	-87.73
+15978.515652	-87.94
+15984.375027	-89.24
+15990.234402	-89.09
+15996.093777	-87.97
+16001.953152	-85.99
+16007.812527	-85.55
+16013.671902	-87.10
+16019.531277	-87.77
+16025.390652	-89.13
+16031.250027	-88.56
+16037.109402	-86.71
+16042.968777	-87.60
+16048.828152	-89.40
+16054.687527	-88.57
+16060.546902	-88.73
+16066.406277	-89.25
+16072.265652	-88.23
+16078.125027	-88.32
+16083.984402	-89.72
+16089.843777	-90.75
+16095.703152	-89.48
+16101.562527	-89.53
+16107.421902	-90.43
+16113.281278	-86.97
+16119.140653	-85.68
+16125.000028	-87.03
+16130.859403	-88.90
+16136.718778	-90.47
+16142.578153	-88.17
+16148.437528	-86.08
+16154.296903	-87.40
+16160.156278	-89.87
+16166.015653	-88.37
+16171.875028	-87.85
+16177.734403	-90.15
+16183.593778	-91.36
+16189.453153	-88.77
+16195.312528	-87.51
+16201.171903	-86.04
+16207.031278	-85.12
+16212.890653	-85.79
+16218.750028	-87.23
+16224.609403	-88.62
+16230.468778	-89.54
+16236.328153	-90.03
+16242.187528	-89.75
+16248.046903	-91.31
+16253.906278	-91.08
+16259.765653	-87.35
+16265.625028	-85.03
+16271.484403	-86.21
+16277.343778	-89.43
+16283.203153	-89.57
+16289.062528	-89.42
+16294.921903	-87.69
+16300.781278	-87.95
+16306.640653	-89.98
+16312.500028	-89.21
+16318.359403	-88.91
+16324.218778	-89.94
+16330.078153	-90.32
+16335.937528	-89.40
+16341.796903	-87.42
+16347.656278	-86.40
+16353.515653	-87.45
+16359.375028	-89.13
+16365.234403	-90.96
+16371.093778	-93.09
+16376.953153	-91.41
+16382.812528	-88.73
+16388.671903	-87.11
+16394.531278	-87.02
+16400.390653	-87.92
+16406.250028	-89.37
+16412.109403	-88.60
+16417.968778	-87.62
+16423.828153	-89.17
+16429.687528	-89.81
+16435.546903	-88.17
+16441.406278	-87.15
+16447.265653	-86.21
+16453.125028	-86.60
+16458.984403	-88.15
+16464.843778	-88.81
+16470.703153	-87.99
+16476.562528	-87.72
+16482.421903	-89.25
+16488.281278	-91.16
+16494.140653	-91.19
+16500.000028	-89.89
+16505.859403	-87.88
+16511.718778	-87.47
+16517.578153	-88.02
+16523.437528	-87.05
+16529.296903	-87.71
+16535.156278	-90.49
+16541.015653	-93.14
+16546.875028	-93.52
+16552.734403	-90.85
+16558.593778	-90.02
+16564.453153	-90.86
+16570.312528	-91.98
+16576.171903	-91.18
+16582.031278	-88.08
+16587.890653	-88.49
+16593.750028	-91.66
+16599.609403	-90.13
+16605.468778	-89.50
+16611.328153	-91.49
+16617.187528	-92.25
+16623.046903	-93.38
+16628.906278	-91.80
+16634.765653	-88.83
+16640.625028	-89.90
+16646.484403	-92.75
+16652.343778	-93.13
+16658.203153	-93.02
+16664.062528	-92.32
+16669.921903	-91.55
+16675.781278	-92.50
+16681.640653	-93.93
+16687.500028	-94.03
+16693.359403	-92.88
+16699.218779	-90.35
+16705.078154	-89.19
+16710.937529	-89.79
+16716.796904	-88.06
+16722.656279	-86.87
+16728.515654	-86.96
+16734.375029	-84.86
+16740.234404	-85.58
+16746.093779	-88.32
+16751.953154	-88.11
+16757.812529	-89.21
+16763.671904	-90.33
+16769.531279	-91.75
+16775.390654	-89.20
+16781.250029	-87.54
+16787.109404	-86.37
+16792.968779	-86.04
+16798.828154	-86.65
+16804.687529	-88.73
+16810.546904	-92.08
+16816.406279	-93.69
+16822.265654	-94.21
+16828.125029	-94.08
+16833.984404	-92.50
+16839.843779	-91.65
+16845.703154	-91.98
+16851.562529	-91.49
+16857.421904	-90.56
+16863.281279	-89.25
+16869.140654	-88.64
+16875.000029	-90.58
+16880.859404	-92.76
+16886.718779	-92.93
+16892.578154	-95.17
+16898.437529	-93.76
+16904.296904	-93.26
+16910.156279	-96.48
+16916.015654	-97.81
+16921.875029	-94.82
+16927.734404	-91.89
+16933.593779	-90.82
+16939.453154	-90.53
+16945.312529	-91.76
+16951.171904	-93.22
+16957.031279	-92.65
+16962.890654	-93.87
+16968.750029	-93.71
+16974.609404	-91.82
+16980.468779	-91.89
+16986.328154	-93.73
+16992.187529	-92.85
+16998.046904	-91.70
+17003.906279	-92.09
+17009.765654	-93.93
+17015.625029	-96.34
+17021.484404	-94.92
+17027.343779	-93.28
+17033.203154	-92.60
+17039.062529	-91.75
+17044.921904	-92.82
+17050.781279	-90.80
+17056.640654	-89.91
+17062.500029	-93.61
+17068.359404	-94.40
+17074.218779	-92.68
+17080.078154	-92.03
+17085.937529	-92.09
+17091.796904	-90.77
+17097.656279	-89.38
+17103.515654	-89.20
+17109.375029	-90.97
+17115.234404	-90.34
+17121.093779	-89.47
+17126.953154	-90.80
+17132.812529	-88.91
+17138.671904	-90.57
+17144.531279	-93.84
+17150.390654	-94.25
+17156.250029	-92.73
+17162.109404	-93.23
+17167.968779	-94.53
+17173.828154	-93.71
+17179.687529	-92.10
+17185.546904	-92.78
+17191.406279	-93.51
+17197.265654	-92.95
+17203.125029	-92.16
+17208.984404	-93.03
+17214.843779	-92.70
+17220.703154	-91.81
+17226.562529	-92.07
+17232.421904	-91.74
+17238.281279	-92.28
+17244.140654	-94.64
+17250.000029	-96.26
+17255.859404	-95.12
+17261.718779	-94.35
+17267.578154	-93.55
+17273.437529	-92.98
+17279.296904	-95.09
+17285.156279	-97.37
+17291.015655	-96.56
+17296.87503	-95.36
+17302.734405	-93.55
+17308.59378	-91.69
+17314.453155	-93.15
+17320.31253	-93.06
+17326.171905	-92.62
+17332.03128	-94.04
+17337.890655	-92.12
+17343.75003	-91.83
+17349.609405	-92.44
+17355.46878	-91.57
+17361.328155	-90.62
+17367.18753	-91.47
+17373.046905	-89.38
+17378.90628	-89.96
+17384.765655	-92.54
+17390.62503	-91.57
+17396.484405	-93.43
+17402.34378	-92.75
+17408.203155	-92.62
+17414.06253	-93.55
+17419.921905	-93.84
+17425.78128	-91.55
+17431.640655	-91.98
+17437.50003	-92.77
+17443.359405	-91.73
+17449.21878	-92.44
+17455.078155	-90.28
+17460.93753	-88.72
+17466.796905	-91.57
+17472.65628	-95.24
+17478.515655	-92.00
+17484.37503	-90.33
+17490.234405	-92.58
+17496.09378	-94.03
+17501.953155	-93.83
+17507.81253	-93.40
+17513.671905	-92.75
+17519.53128	-91.16
+17525.390655	-90.45
+17531.25003	-92.06
+17537.109405	-94.19
+17542.96878	-91.83
+17548.828155	-90.73
+17554.68753	-91.76
+17560.546905	-91.27
+17566.40628	-89.88
+17572.265655	-90.44
+17578.12503	-90.74
+17583.984405	-90.51
+17589.84378	-90.66
+17595.703155	-90.75
+17601.56253	-92.15
+17607.421905	-92.26
+17613.28128	-92.40
+17619.140655	-92.99
+17625.00003	-92.26
+17630.859405	-91.43
+17636.71878	-90.34
+17642.578155	-88.13
+17648.43753	-87.84
+17654.296905	-90.39
+17660.15628	-93.24
+17666.015655	-91.41
+17671.87503	-90.96
+17677.734405	-91.52
+17683.59378	-94.00
+17689.453155	-92.16
+17695.31253	-89.55
+17701.171905	-90.56
+17707.03128	-93.21
+17712.890655	-90.84
+17718.75003	-90.10
+17724.609405	-92.59
+17730.46878	-90.08
+17736.328155	-87.21
+17742.18753	-87.84
+17748.046905	-91.22
+17753.90628	-90.57
+17759.765655	-89.14
+17765.62503	-90.12
+17771.484405	-91.94
+17777.34378	-91.21
+17783.203155	-92.17
+17789.06253	-93.31
+17794.921905	-92.75
+17800.78128	-92.71
+17806.640655	-91.75
+17812.50003	-90.53
+17818.359405	-91.39
+17824.21878	-92.10
+17830.078155	-90.87
+17835.93753	-91.20
+17841.796905	-90.75
+17847.65628	-91.83
+17853.515655	-89.98
+17859.37503	-88.01
+17865.234405	-90.83
+17871.09378	-89.80
+17876.953156	-88.34
+17882.812531	-91.37
+17888.671906	-93.85
+17894.531281	-93.52
+17900.390656	-94.18
+17906.250031	-93.85
+17912.109406	-91.09
+17917.968781	-92.20
+17923.828156	-94.39
+17929.687531	-92.59
+17935.546906	-90.62
+17941.406281	-90.87
+17947.265656	-93.12
+17953.125031	-92.14
+17958.984406	-88.41
+17964.843781	-89.72
+17970.703156	-94.67
+17976.562531	-96.65
+17982.421906	-95.52
+17988.281281	-94.11
+17994.140656	-94.97
+18000.000031	-95.78
+18005.859406	-95.48
+18011.718781	-95.11
+18017.578156	-95.17
+18023.437531	-97.96
+18029.296906	-97.37
+18035.156281	-94.86
+18041.015656	-94.82
+18046.875031	-96.89
+18052.734406	-95.97
+18058.593781	-94.73
+18064.453156	-92.90
+18070.312531	-90.59
+18076.171906	-93.04
+18082.031281	-97.04
+18087.890656	-96.04
+18093.750031	-96.27
+18099.609406	-94.65
+18105.468781	-91.63
+18111.328156	-92.43
+18117.187531	-93.33
+18123.046906	-93.58
+18128.906281	-96.14
+18134.765656	-93.77
+18140.625031	-92.12
+18146.484406	-93.03
+18152.343781	-94.16
+18158.203156	-96.49
+18164.062531	-98.64
+18169.921906	-96.96
+18175.781281	-97.89
+18181.640656	-98.71
+18187.500031	-97.57
+18193.359406	-97.69
+18199.218781	-97.83
+18205.078156	-96.75
+18210.937531	-97.28
+18216.796906	-96.94
+18222.656281	-96.47
+18228.515656	-95.62
+18234.375031	-94.79
+18240.234406	-93.82
+18246.093781	-92.71
+18251.953156	-93.78
+18257.812531	-95.20
+18263.671906	-96.54
+18269.531281	-97.25
+18275.390656	-93.46
+18281.250031	-92.33
+18287.109406	-94.61
+18292.968781	-95.64
+18298.828156	-95.34
+18304.687531	-97.81
+18310.546906	-98.64
+18316.406281	-97.22
+18322.265656	-96.06
+18328.125031	-93.76
+18333.984406	-91.70
+18339.843781	-92.22
+18345.703156	-94.49
+18351.562531	-94.43
+18357.421906	-93.55
+18363.281281	-95.58
+18369.140656	-98.00
+18375.000031	-96.62
+18380.859406	-95.75
+18386.718781	-96.02
+18392.578156	-97.98
+18398.437531	-98.69
+18404.296906	-96.50
+18410.156281	-96.46
+18416.015656	-97.50
+18421.875031	-99.74
+18427.734406	-98.72
+18433.593781	-97.89
+18439.453156	-96.86
+18445.312531	-95.84
+18451.171906	-94.56
+18457.031282	-94.35
+18462.890657	-95.49
+18468.750032	-96.11
+18474.609407	-95.11
+18480.468782	-92.79
+18486.328157	-92.87
+18492.187532	-95.41
+18498.046907	-97.08
+18503.906282	-98.73
+18509.765657	-99.18
+18515.625032	-97.08
+18521.484407	-95.38
+18527.343782	-94.22
+18533.203157	-95.57
+18539.062532	-94.92
+18544.921907	-95.90
+18550.781282	-94.93
+18556.640657	-95.29
+18562.500032	-94.60
+18568.359407	-93.30
+18574.218782	-94.36
+18580.078157	-93.17
+18585.937532	-93.98
+18591.796907	-96.16
+18597.656282	-95.89
+18603.515657	-96.40
+18609.375032	-97.47
+18615.234407	-94.78
+18621.093782	-92.89
+18626.953157	-94.14
+18632.812532	-96.83
+18638.671907	-96.45
+18644.531282	-94.76
+18650.390657	-91.75
+18656.250032	-91.92
+18662.109407	-95.96
+18667.968782	-99.31
+18673.828157	-99.66
+18679.687532	-101.54
+18685.546907	-98.49
+18691.406282	-97.02
+18697.265657	-99.53
+18703.125032	-100.11
+18708.984407	-100.45
+18714.843782	-99.12
+18720.703157	-93.81
+18726.562532	-92.92
+18732.421907	-96.55
+18738.281282	-96.92
+18744.140657	-95.18
+18750.000032	-96.29
+18755.859407	-98.32
+18761.718782	-97.76
+18767.578157	-97.23
+18773.437532	-97.21
+18779.296907	-97.57
+18785.156282	-97.43
+18791.015657	-95.88
+18796.875032	-94.85
+18802.734407	-95.77
+18808.593782	-97.45
+18814.453157	-96.54
+18820.312532	-96.26
+18826.171907	-95.28
+18832.031282	-92.22
+18837.890657	-91.82
+18843.750032	-93.54
+18849.609407	-94.27
+18855.468782	-94.62
+18861.328157	-95.18
+18867.187532	-96.60
+18873.046907	-95.05
+18878.906282	-93.48
+18884.765657	-96.48
+18890.625032	-97.24
+18896.484407	-97.58
+18902.343782	-97.40
+18908.203157	-95.87
+18914.062532	-95.02
+18919.921907	-94.99
+18925.781282	-96.00
+18931.640657	-96.18
+18937.500032	-92.68
+18943.359407	-91.34
+18949.218782	-94.72
+18955.078157	-97.81
+18960.937532	-95.73
+18966.796907	-95.59
+18972.656282	-93.44
+18978.515657	-92.18
+18984.375032	-93.97
+18990.234407	-96.51
+18996.093782	-95.36
+19001.953157	-93.51
+19007.812532	-93.80
+19013.671907	-92.82
+19019.531282	-92.50
+19025.390657	-95.35
+19031.250032	-95.80
+19037.109407	-93.38
+19042.968783	-94.08
+19048.828158	-95.91
+19054.687533	-95.20
+19060.546908	-95.07
+19066.406283	-94.98
+19072.265658	-93.84
+19078.125033	-94.77
+19083.984408	-96.53
+19089.843783	-93.29
+19095.703158	-92.33
+19101.562533	-94.85
+19107.421908	-97.05
+19113.281283	-96.76
+19119.140658	-95.52
+19125.000033	-96.13
+19130.859408	-97.08
+19136.718783	-96.37
+19142.578158	-96.10
+19148.437533	-95.73
+19154.296908	-93.97
+19160.156283	-94.39
+19166.015658	-94.20
+19171.875033	-93.43
+19177.734408	-94.39
+19183.593783	-94.70
+19189.453158	-94.71
+19195.312533	-94.59
+19201.171908	-94.42
+19207.031283	-93.96
+19212.890658	-94.83
+19218.750033	-97.38
+19224.609408	-95.39
+19230.468783	-93.16
+19236.328158	-95.32
+19242.187533	-94.62
+19248.046908	-90.65
+19253.906283	-90.13
+19259.765658	-93.09
+19265.625033	-96.16
+19271.484408	-97.76
+19277.343783	-98.50
+19283.203158	-96.31
+19289.062533	-94.58
+19294.921908	-95.52
+19300.781283	-92.99
+19306.640658	-92.30
+19312.500033	-95.44
+19318.359408	-97.61
+19324.218783	-95.75
+19330.078158	-93.57
+19335.937533	-91.45
+19341.796908	-91.34
+19347.656283	-93.32
+19353.515658	-96.41
+19359.375033	-97.80
+19365.234408	-97.06
+19371.093783	-95.69
+19376.953158	-96.01
+19382.812533	-95.19
+19388.671908	-93.49
+19394.531283	-93.64
+19400.390658	-94.52
+19406.250033	-96.51
+19412.109408	-98.12
+19417.968783	-95.86
+19423.828158	-94.21
+19429.687533	-96.05
+19435.546908	-97.37
+19441.406283	-96.15
+19447.265658	-95.73
+19453.125033	-94.45
+19458.984408	-95.52
+19464.843783	-97.07
+19470.703158	-94.37
+19476.562533	-95.45
+19482.421908	-96.98
+19488.281283	-94.79
+19494.140658	-95.52
+19500.000033	-98.89
+19505.859408	-101.21
+19511.718783	-98.51
+19517.578158	-97.45
+19523.437533	-95.24
+19529.296908	-93.47
+19535.156283	-96.28
+19541.015658	-99.37
+19546.875033	-95.04
+19552.734408	-93.98
+19558.593783	-96.34
+19564.453158	-98.61
+19570.312533	-98.28
+19576.171908	-96.80
+19582.031283	-95.48
+19587.890658	-96.47
+19593.750033	-98.00
+19599.609408	-95.68
+19605.468783	-96.98
+19611.328158	-97.96
+19617.187533	-96.66
+19623.046908	-96.10
+19628.906284	-94.65
+19634.765659	-94.57
+19640.625034	-95.62
+19646.484409	-93.85
+19652.343784	-92.84
+19658.203159	-94.37
+19664.062534	-97.53
+19669.921909	-98.57
+19675.781284	-97.70
+19681.640659	-95.93
+19687.500034	-96.85
+19693.359409	-93.54
+19699.218784	-92.18
+19705.078159	-95.34
+19710.937534	-95.48
+19716.796909	-94.14
+19722.656284	-96.69
+19728.515659	-99.83
+19734.375034	-98.95
+19740.234409	-96.50
+19746.093784	-95.51
+19751.953159	-96.42
+19757.812534	-97.70
+19763.671909	-97.90
+19769.531284	-99.74
+19775.390659	-100.42
+19781.250034	-99.34
+19787.109409	-96.47
+19792.968784	-95.91
+19798.828159	-96.74
+19804.687534	-97.29
+19810.546909	-96.92
+19816.406284	-96.65
+19822.265659	-97.66
+19828.125034	-98.39
+19833.984409	-98.86
+19839.843784	-99.22
+19845.703159	-95.65
+19851.562534	-94.19
+19857.421909	-96.98
+19863.281284	-100.88
+19869.140659	-100.56
+19875.000034	-99.15
+19880.859409	-96.24
+19886.718784	-93.01
+19892.578159	-92.10
+19898.437534	-94.21
+19904.296909	-95.82
+19910.156284	-95.11
+19916.015659	-93.80
+19921.875034	-91.04
+19927.734409	-92.75
+19933.593784	-95.31
+19939.453159	-94.65
+19945.312534	-97.12
+19951.171909	-98.63
+19957.031284	-96.88
+19962.890659	-96.76
+19968.750034	-98.25
+19974.609409	-98.80
+19980.468784	-98.17
+19986.328159	-96.66
+19992.187534	-94.76
+19998.046909	-95.06
+20003.906284	-97.13
+20009.765659	-97.76
+20015.625034	-96.49
+20021.484409	-96.31
+20027.343784	-96.99
+20033.203159	-98.04
+20039.062534	-98.30
+20044.921909	-98.72
+20050.781284	-94.12
+20056.640659	-92.84
+20062.500034	-95.49
+20068.359409	-96.25
+20074.218784	-98.66
+20080.078159	-98.48
+20085.937534	-95.65
+20091.796909	-94.12
+20097.656284	-95.03
+20103.515659	-97.97
+20109.375034	-97.16
+20115.234409	-93.70
+20121.093784	-93.04
+20126.953159	-94.29
+20132.812534	-96.05
+20138.671909	-95.44
+20144.531284	-93.41
+20150.390659	-95.03
+20156.250034	-93.92
+20162.109409	-90.37
+20167.968784	-90.29
+20173.828159	-92.34
+20179.687534	-95.78
+20185.546909	-95.88
+20191.406284	-95.20
+20197.265659	-96.79
+20203.125034	-97.22
+20208.984409	-97.08
+20214.843785	-99.10
+20220.70316	-99.79
+20226.562535	-98.82
+20232.42191	-96.41
+20238.281285	-94.64
+20244.14066	-91.39
+20250.000035	-90.34
+20255.85941	-94.03
+20261.718785	-95.55
+20267.57816	-95.48
+20273.437535	-95.46
+20279.29691	-94.55
+20285.156285	-94.82
+20291.01566	-91.44
+20296.875035	-91.79
+20302.73441	-94.31
+20308.593785	-91.42
+20314.45316	-89.80
+20320.312535	-91.72
+20326.17191	-96.09
+20332.031285	-97.12
+20337.89066	-96.93
+20343.750035	-94.07
+20349.60941	-91.95
+20355.468785	-93.75
+20361.32816	-96.61
+20367.187535	-93.29
+20373.04691	-93.29
+20378.906285	-96.74
+20384.76566	-96.91
+20390.625035	-97.24
+20396.48441	-98.74
+20402.343785	-97.87
+20408.20316	-95.23
+20414.062535	-93.29
+20419.92191	-92.83
+20425.781285	-93.92
+20431.64066	-96.67
+20437.500035	-100.11
+20443.35941	-98.64
+20449.218785	-96.57
+20455.07816	-95.74
+20460.937535	-95.96
+20466.79691	-97.01
+20472.656285	-97.16
+20478.51566	-95.21
+20484.375035	-95.20
+20490.23441	-99.12
+20496.093785	-98.70
+20501.95316	-94.92
+20507.812535	-94.56
+20513.67191	-97.12
+20519.531285	-98.22
+20525.39066	-95.40
+20531.250035	-95.88
+20537.10941	-98.23
+20542.968785	-99.14
+20548.82816	-100.58
+20554.687535	-101.82
+20560.54691	-104.86
+20566.406285	-105.57
+20572.26566	-99.84
+20578.125035	-96.74
+20583.98441	-98.14
+20589.843785	-96.37
+20595.70316	-93.90
+20601.562535	-94.91
+20607.42191	-96.29
+20613.281285	-96.29
+20619.14066	-98.25
+20625.000035	-99.32
+20630.85941	-99.04
+20636.718785	-100.76
+20642.57816	-98.68
+20648.437535	-98.03
+20654.29691	-99.44
+20660.156285	-102.02
+20666.01566	-100.40
+20671.875035	-100.34
+20677.73441	-100.96
+20683.593785	-99.27
+20689.45316	-99.01
+20695.312535	-101.30
+20701.17191	-98.87
+20707.031285	-94.80
+20712.89066	-94.17
+20718.750035	-96.44
+20724.60941	-97.96
+20730.468785	-100.74
+20736.32816	-100.79
+20742.187535	-99.03
+20748.04691	-96.88
+20753.906285	-97.10
+20759.76566	-98.79
+20765.625035	-95.79
+20771.48441	-96.14
+20777.343785	-98.58
+20783.20316	-97.26
+20789.062535	-96.33
+20794.92191	-96.44
+20800.781286	-97.92
+20806.640661	-100.68
+20812.500036	-102.65
+20818.359411	-103.20
+20824.218786	-100.20
+20830.078161	-99.37
+20835.937536	-99.77
+20841.796911	-97.17
+20847.656286	-97.01
+20853.515661	-93.28
+20859.375036	-92.09
+20865.234411	-95.18
+20871.093786	-98.73
+20876.953161	-99.99
+20882.812536	-102.01
+20888.671911	-99.61
+20894.531286	-99.48
+20900.390661	-102.43
+20906.250036	-101.60
+20912.109411	-100.33
+20917.968786	-101.86
+20923.828161	-103.30
+20929.687536	-103.76
+20935.546911	-101.87
+20941.406286	-100.46
+20947.265661	-97.82
+20953.125036	-95.86
+20958.984411	-95.74
+20964.843786	-94.91
+20970.703161	-96.75
+20976.562536	-101.37
+20982.421911	-103.98
+20988.281286	-104.89
+20994.140661	-101.45
+21000.000036	-97.25
+21005.859411	-94.11
+21011.718786	-93.90
+21017.578161	-95.47
+21023.437536	-96.97
+21029.296911	-101.12
+21035.156286	-102.45
+21041.015661	-102.94
+21046.875036	-103.75
+21052.734411	-100.65
+21058.593786	-99.01
+21064.453161	-96.58
+21070.312536	-95.92
+21076.171911	-95.87
+21082.031286	-95.79
+21087.890661	-97.86
+21093.750036	-100.10
+21099.609411	-100.13
+21105.468786	-98.05
+21111.328161	-97.40
+21117.187536	-97.71
+21123.046911	-97.81
+21128.906286	-99.58
+21134.765661	-103.56
+21140.625036	-102.57
+21146.484411	-101.59
+21152.343786	-100.90
+21158.203161	-98.08
+21164.062536	-98.06
+21169.921911	-101.31
+21175.781286	-101.33
+21181.640661	-102.44
+21187.500036	-100.38
+21193.359411	-95.19
+21199.218786	-93.49
+21205.078161	-95.69
+21210.937536	-98.24
+21216.796911	-95.44
+21222.656286	-95.60
+21228.515661	-98.03
+21234.375036	-96.62
+21240.234411	-96.67
+21246.093786	-98.92
+21251.953161	-97.85
+21257.812536	-99.73
+21263.671911	-100.64
+21269.531286	-98.87
+21275.390661	-100.77
+21281.250036	-98.59
+21287.109411	-98.13
+21292.968786	-100.99
+21298.828161	-100.80
+21304.687536	-102.62
+21310.546911	-102.80
+21316.406286	-102.48
+21322.265661	-102.00
+21328.125036	-101.70
+21333.984411	-101.60
+21339.843786	-101.89
+21345.703161	-100.20
+21351.562536	-98.07
+21357.421911	-97.53
+21363.281286	-100.12
+21369.140661	-99.27
+21375.000036	-99.52
+21380.859411	-101.06
+21386.718787	-98.63
+21392.578162	-98.30
+21398.437537	-97.16
+21404.296912	-97.07
+21410.156287	-99.05
+21416.015662	-100.48
+21421.875037	-100.87
+21427.734412	-102.41
+21433.593787	-104.07
+21439.453162	-101.85
+21445.312537	-101.72
+21451.171912	-103.80
+21457.031287	-102.82
+21462.890662	-101.90
+21468.750037	-99.88
+21474.609412	-97.52
+21480.468787	-96.89
+21486.328162	-97.23
+21492.187537	-99.02
+21498.046912	-100.25
+21503.906287	-97.76
+21509.765662	-98.45
+21515.625037	-100.63
+21521.484412	-100.00
+21527.343787	-100.05
+21533.203162	-101.31
+21539.062537	-101.06
+21544.921912	-98.59
+21550.781287	-98.39
+21556.640662	-99.10
+21562.500037	-96.26
+21568.359412	-94.05
+21574.218787	-94.63
+21580.078162	-97.14
+21585.937537	-101.03
+21591.796912	-102.87
+21597.656287	-100.40
+21603.515662	-99.57
+21609.375037	-100.82
+21615.234412	-102.38
+21621.093787	-104.96
+21626.953162	-102.91
+21632.812537	-101.09
+21638.671912	-102.51
+21644.531287	-100.76
+21650.390662	-97.93
+21656.250037	-99.54
+21662.109412	-103.35
+21667.968787	-101.86
+21673.828162	-98.63
+21679.687537	-98.38
+21685.546912	-101.13
+21691.406287	-103.78
+21697.265662	-103.39
+21703.125037	-103.20
+21708.984412	-101.51
+21714.843787	-99.82
+21720.703162	-101.94
+21726.562537	-103.71
+21732.421912	-101.43
+21738.281287	-101.65
+21744.140662	-103.00
+21750.000037	-101.35
+21755.859412	-99.92
+21761.718787	-100.53
+21767.578162	-100.37
+21773.437537	-101.03
+21779.296912	-102.79
+21785.156287	-100.82
+21791.015662	-100.50
+21796.875037	-103.10
+21802.734412	-103.61
+21808.593787	-102.27
+21814.453162	-103.58
+21820.312537	-104.11
+21826.171912	-103.35
+21832.031287	-102.62
+21837.890662	-101.96
+21843.750037	-101.43
+21849.609412	-100.11
+21855.468787	-100.43
+21861.328162	-102.80
+21867.187537	-101.39
+21873.046912	-100.14
+21878.906287	-102.03
+21884.765662	-103.05
+21890.625037	-101.36
+21896.484412	-102.78
+21902.343787	-100.98
+21908.203162	-99.35
+21914.062537	-100.96
+21919.921912	-100.65
+21925.781287	-101.29
+21931.640662	-102.08
+21937.500037	-99.40
+21943.359412	-100.69
+21949.218787	-104.68
+21955.078162	-102.72
+21960.937537	-103.38
+21966.796912	-104.89
+21972.656288	-104.48
+21978.515663	-103.27
+21984.375038	-98.92
+21990.234413	-97.54
+21996.093788	-98.08
+22001.953163	-96.08
+22007.812538	-97.64
+22013.671913	-101.50
+22019.531288	-101.92
+22025.390663	-101.51
+22031.250038	-99.36
+22037.109413	-99.09
+22042.968788	-101.02
+22048.828163	-103.45
+22054.687538	-100.18
+22060.546913	-98.78
+22066.406288	-100.10
+22072.265663	-99.68
+22078.125038	-97.54
+22083.984413	-98.03
+22089.843788	-101.50
+22095.703163	-102.80
+22101.562538	-103.30
+22107.421913	-105.79
+22113.281288	-107.06
+22119.140663	-105.54
+22125.000038	-101.31
+22130.859413	-98.72
+22136.718788	-100.67
+22142.578163	-100.02
+22148.437538	-98.79
+22154.296913	-101.22
+22160.156288	-101.59
+22166.015663	-101.69
+22171.875038	-102.99
+22177.734413	-104.37
+22183.593788	-104.53
+22189.453163	-102.48
+22195.312538	-102.08
+22201.171913	-103.03
+22207.031288	-102.19
+22212.890663	-102.17
+22218.750038	-103.88
+22224.609413	-106.11
+22230.468788	-106.44
+22236.328163	-106.17
+22242.187538	-104.70
+22248.046913	-105.63
+22253.906288	-105.82
+22259.765663	-102.71
+22265.625038	-99.98
+22271.484413	-99.57
+22277.343788	-102.20
+22283.203163	-104.79
+22289.062538	-106.38
+22294.921913	-104.11
+22300.781288	-101.49
+22306.640663	-102.87
+22312.500038	-104.14
+22318.359413	-100.45
+22324.218788	-100.16
+22330.078163	-101.86
+22335.937538	-102.57
+22341.796913	-105.43
+22347.656288	-106.39
+22353.515663	-103.58
+22359.375038	-102.03
+22365.234413	-101.99
+22371.093788	-101.92
+22376.953163	-102.22
+22382.812538	-103.72
+22388.671913	-103.51
+22394.531288	-102.29
+22400.390663	-102.74
+22406.250038	-105.89
+22412.109413	-109.09
+22417.968788	-107.29
+22423.828163	-103.65
+22429.687538	-102.76
+22435.546913	-104.59
+22441.406288	-105.19
+22447.265663	-104.71
+22453.125038	-104.59
+22458.984413	-105.07
+22464.843788	-101.76
+22470.703163	-100.29
+22476.562538	-101.85
+22482.421913	-104.90
+22488.281288	-106.12
+22494.140663	-105.60
+22500.000038	-105.07
+22505.859413	-106.02
+22511.718788	-106.36
+22517.578163	-107.55
+22523.437538	-108.74
+22529.296913	-107.60
+22535.156288	-106.84
+22541.015663	-101.74
+22546.875038	-100.27
+22552.734413	-103.37
+22558.593789	-106.29
+22564.453164	-106.45
+22570.312539	-106.72
+22576.171914	-107.63
+22582.031289	-107.14
+22587.890664	-107.01
+22593.750039	-106.41
+22599.609414	-105.84
+22605.468789	-107.14
+22611.328164	-106.76
+22617.187539	-107.12
+22623.046914	-105.79
+22628.906289	-105.11
+22634.765664	-105.88
+22640.625039	-106.40
+22646.484414	-107.59
+22652.343789	-107.44
+22658.203164	-106.01
+22664.062539	-104.08
+22669.921914	-102.49
+22675.781289	-102.81
+22681.640664	-105.28
+22687.500039	-104.70
+22693.359414	-103.24
+22699.218789	-105.22
+22705.078164	-107.68
+22710.937539	-107.11
+22716.796914	-105.59
+22722.656289	-103.19
+22728.515664	-102.25
+22734.375039	-104.17
+22740.234414	-105.26
+22746.093789	-104.62
+22751.953164	-107.35
+22757.812539	-108.56
+22763.671914	-104.14
+22769.531289	-101.97
+22775.390664	-102.36
+22781.250039	-101.96
+22787.109414	-103.66
+22792.968789	-104.82
+22798.828164	-102.74
+22804.687539	-103.45
+22810.546914	-104.75
+22816.406289	-104.56
+22822.265664	-105.07
+22828.125039	-107.02
+22833.984414	-102.52
+22839.843789	-98.53
+22845.703164	-99.22
+22851.562539	-103.20
+22857.421914	-103.17
+22863.281289	-104.87
+22869.140664	-104.59
+22875.000039	-102.65
+22880.859414	-103.92
+22886.718789	-107.45
+22892.578164	-109.99
+22898.437539	-108.70
+22904.296914	-107.57
+22910.156289	-107.11
+22916.015664	-105.53
+22921.875039	-104.67
+22927.734414	-106.33
+22933.593789	-107.18
+22939.453164	-105.38
+22945.312539	-105.42
+22951.171914	-105.54
+22957.031289	-105.14
+22962.890664	-103.90
+22968.750039	-104.86
+22974.609414	-107.20
+22980.468789	-106.66
+22986.328164	-104.73
+22992.187539	-102.00
+22998.046914	-100.84
+23003.906289	-101.57
+23009.765664	-103.88
+23015.625039	-106.79
+23021.484414	-106.88
+23027.343789	-104.88
+23033.203164	-102.76
+23039.062539	-103.23
+23044.921914	-106.35
+23050.781289	-106.54
+23056.640664	-106.25
+23062.500039	-106.20
+23068.359414	-105.95
+23074.218789	-105.38
+23080.078164	-105.50
+23085.937539	-106.33
+23091.796914	-105.71
+23097.656289	-103.74
+23103.515664	-103.70
+23109.375039	-104.93
+23115.234414	-105.76
+23121.093789	-107.90
+23126.953164	-110.32
+23132.812539	-107.52
+23138.671914	-102.98
+23144.531289	-103.33
+23150.390665	-106.21
+23156.25004	-105.18
+23162.109415	-106.08
+23167.96879	-106.85
+23173.828165	-104.18
+23179.68754	-102.51
+23185.546915	-103.45
+23191.40629	-106.38
+23197.265665	-107.58
+23203.12504	-106.72
+23208.984415	-105.98
+23214.84379	-106.85
+23220.703165	-107.54
+23226.56254	-104.99
+23232.421915	-102.33
+23238.28129	-100.97
+23244.140665	-101.21
+23250.00004	-103.31
+23255.859415	-106.23
+23261.71879	-109.10
+23267.578165	-110.58
+23273.43754	-110.77
+23279.296915	-110.39
+23285.15629	-108.73
+23291.015665	-106.98
+23296.87504	-105.74
+23302.734415	-104.34
+23308.59379	-105.14
+23314.453165	-106.64
+23320.31254	-105.77
+23326.171915	-106.23
+23332.03129	-107.41
+23337.890665	-105.61
+23343.75004	-105.66
+23349.609415	-107.90
+23355.46879	-109.38
+23361.328165	-109.77
+23367.18754	-109.05
+23373.046915	-109.01
+23378.90629	-107.53
+23384.765665	-106.08
+23390.62504	-106.94
+23396.484415	-107.47
+23402.34379	-107.51
+23408.203165	-109.29
+23414.06254	-108.84
+23419.921915	-107.16
+23425.78129	-105.95
+23431.640665	-106.35
+23437.50004	-107.48
+23443.359415	-106.85
+23449.21879	-107.24
+23455.078165	-105.31
+23460.93754	-102.99
+23466.796915	-105.02
+23472.65629	-109.20
+23478.515665	-109.90
+23484.37504	-110.34
+23490.234415	-109.90
+23496.09379	-107.63
+23501.953165	-107.01
+23507.81254	-108.81
+23513.671915	-108.46
+23519.53129	-108.85
+23525.390665	-109.29
+23531.25004	-107.57
+23537.109415	-108.03
+23542.96879	-110.33
+23548.828165	-110.14
+23554.68754	-110.25
+23560.546915	-108.32
+23566.40629	-105.99
+23572.265665	-107.30
+23578.12504	-109.97
+23583.984415	-106.99
+23589.84379	-105.38
+23595.703165	-106.30
+23601.56254	-107.59
+23607.421915	-106.23
+23613.28129	-104.22
+23619.140665	-104.54
+23625.00004	-105.75
+23630.859415	-107.72
+23636.71879	-110.00
+23642.578165	-110.97
+23648.43754	-111.86
+23654.296915	-111.14
+23660.15629	-110.11
+23666.015665	-109.23
+23671.87504	-109.53
+23677.734415	-108.88
+23683.59379	-108.16
+23689.453165	-109.64
+23695.31254	-108.62
+23701.171915	-107.40
+23707.03129	-105.41
+23712.890665	-102.23
+23718.75004	-102.26
+23724.609415	-105.93
+23730.46879	-105.89
+23736.328166	-103.73
+23742.187541	-104.88
+23748.046916	-108.43
+23753.906291	-109.06
+23759.765666	-108.47
+23765.625041	-108.76
+23771.484416	-110.95
+23777.343791	-112.03
+23783.203166	-111.22
+23789.062541	-109.51
+23794.921916	-108.15
+23800.781291	-108.75
+23806.640666	-107.92
+23812.500041	-106.13
+23818.359416	-108.50
+23824.218791	-107.21
+23830.078166	-106.42
+23835.937541	-108.29
+23841.796916	-108.24
+23847.656291	-106.35
+23853.515666	-105.74
+23859.375041	-107.54
+23865.234416	-109.70
+23871.093791	-108.39
+23876.953166	-108.40
+23882.812541	-107.44
+23888.671916	-105.52
+23894.531291	-105.89
+23900.390666	-108.77
+23906.250041	-111.48
+23912.109416	-110.25
+23917.968791	-106.68
+23923.828166	-106.31
+23929.687541	-108.91
+23935.546916	-106.48
+23941.406291	-105.46
+23947.265666	-105.82
+23953.125041	-106.95
+23958.984416	-107.84
+23964.843791	-106.15
+23970.703166	-102.82
+23976.562541	-100.04
+23982.421916	-100.23
+23988.281291	-103.27
+23994.140666	-101.77