#!/usr/bin/env python3
# -*- coding: utf-8 -*-
from __future__ import annotations
"""!
Author: J.A. de Jong - ASCEE

Description: Measurement class

The ASCEE hdf5 measurement file format contains the following fields:

- Attributes:

'LASP_VERSION_MAJOR': int The major version of LASP which which the recording has been performed.
'LASP_VERSION_MINOR': int The minor version of LASP which which the recording has been performed.

'samplerate': The audio data sample rate in Hz [float]
'nchannels': The number of audio channels in the file List[float]
'sensitivity': (Optionally) the stored sensitivity of the record channels.
               This can be a single value, or a list of sensitivities for
               each channel. Both representations are allowed. List[float]

For measurement files of LASP < v1.0
'qtys' : (Optionally): list of quantities that is recorded for each channel',
if this array is not found. Quantities are defaulted to 'Number / Full scale'

'type_int': A specified measurement type that can be used programmatically. It can be read out as an enumeration variant of type "MeasurementType". See code below of implemented measurement types.

For measurement files of LASP >= 1.0

- Datasets:

'audio': 3-dimensional array of blocks of audio data. The first axis is the
block index, the second axis the sample number and the third axis is the
channel number. The data type is either int16, int32 or float64 / float32. If
raw data is stored as integer values (int16, int32), the actual values should
be pre-scaled by its maximum positive number (2**(nb-1) - 1), such that the
corresponding 'number' lies between -1.0 and 1.0. To stay backwards-compatible, the dataset is always called 'audio' despite it being possible that other types of data is stored in the dataset (such as voltages, accelerations etc).

'video': 4-dimensional array of video frames. The first index is the frame
         number, the second the x-value of the pixel and the third is the
         y-value of the pixel. Then, the last axis is the color. This axis has
         length 3 and the colors are stored as (r,g,b). Where typically a
         color depth of 256 is used (np.uint8 data format)

The video dataset can possibly be not present in the data.

"""

__all__ = ["Measurement", "scaleBlockSens", "MeasurementType"]
from contextlib import contextmanager
from weakref import WeakValueDictionary
import h5py as h5
import uuid
import pathlib
import glob
import itertools
import numpy as np
from enum import Enum, unique
from .lasp_config import LASP_NUMPY_FLOAT_TYPE
from scipy.io import wavfile
import os, time, wave, logging
from .lasp_common import SIQtys, Qty, getFreq
from .lasp_version import LASP_VERSION_MAJOR, LASP_VERSION_MINOR
from .lasp_cpp import Window, DaqChannel, AvPowerSpectra
from typing import List
from functools import lru_cache

# Measurement file extension
MEXT = 'h5'
DOTMEXT = f'.{MEXT}'


@unique
class MeasurementType(Enum):
    """
    Measurement flags related to the measurement. Stored as bit flags in the measurement file. This is for possible changes in the API later.
    """

    # Not specific measurement type    
    NotSpecific = 0
    
    # Measurement serves as an insertion loss reference measurement
    ILReference = 1 << 0

    # Measurement is general calibration measurement (to calibrate sensor in a certain way)
    CALGeneral = 1 << 1

    # Measurement serves as impedance tube calibration (short tube case / ref. plane at origin)
    muZCalShort = 1 << 2

    # Measurement serves as impedance tube calibration (long tube case)
    muZCalLong = 1 << 3

    # Series impedance reference
    muZSeriesImpedanceRef = 1 << 4

    def __str__(self):
        match self:
            case MeasurementType.NotSpecific: return '-'
            case MeasurementType.ILReference: return 'Insertion loss reference'
            case MeasurementType.CALGeneral: return 'General calibration'
            case MeasurementType.muZCalShort: return 'ASCEE μZ short length calibration'
            case MeasurementType.muZCalLong: return 'ASCEE μZ long length calibration'
            case MeasurementType.muZSeriesImpedanceRef: return 'ASCEE μZ series impedance empty reference'
            case _: raise ValueError("Not a MeasurementType")


def getSampWidth(dtype):
    """Returns the width of a single sample in **bytes**.

    Args:
        dtype: numpy dtype

    Returns:
        Size of a sample in bytes (int)
    """
    if dtype in (np.int32, np.float32):
        return 4
    elif dtype == np.int16:
        return 2
    elif dtype == np.float64:
        return 8
    else:
        raise ValueError("Invalid data type: %s" % dtype)


def scaleBlockSens(block, sens):
    """Scale a block of raw data to return raw acoustic pressure data.

    Args:
        block: block of raw data with integer data type
        sens: array of sensitivity coeficients for
        each channel.
    """
    sens = np.asarray(sens)
    assert sens.size == block.shape[1]
    if np.issubdtype(block.dtype.type, np.integer):
        sw = getSampWidth(block.dtype)
        fac = 2 ** (8 * sw - 1) - 1
    else:
        fac = 1.0
    return block.astype(LASP_NUMPY_FLOAT_TYPE) / fac / sens[np.newaxis,:]


class IterRawData:
    """Iterate over stored blocks if the raw measurement data of a h5 file."""

    def __init__(self, f, channels, **kwargs):
        """Initialize a BlockIter object.

        Args:
            f: Audio dataset in the h5 file, accessed as f['audio']
            channels: list of channel indices to use
            istart: index of first sample
            istop: index of last sample (not including istop)

        """
        assert isinstance(channels, list)
        fa = f["audio"]
        self.fa = fa
        self.i = 0

        nblocks = fa.shape[0]
        blocksize = fa.shape[1]
        self.blocksize = blocksize
        # nchannels = fa.shape[2]
        self.channels = channels

        self.istart = kwargs.pop("istart", 0)
        self.istop = kwargs.pop("istop", blocksize * nblocks)

        self.firstblock = self.istart // blocksize
        self.lastblock = self.istop // blocksize
        if self.istop % blocksize == 0:
            self.lastblock -= 1

        self.firstblock_start_offset = self.istart % blocksize

        if self.istop < 0:
            self.istop += blocksize * nblocks
        if self.istop <= self.istart:
            raise ValueError("Stop index is smaller than start index")

        if self.istop != blocksize * nblocks:
            self.lastblock_stop_offset = self.istop % blocksize
        else:
            self.lastblock_stop_offset = blocksize

    def __iter__(self):
        return self

    def __next__(self):
        """Return the next block."""
        fa = self.fa

        # nblocks_to_return = self.lastblock-self.firstblock+1

        block = self.firstblock + self.i

        if block > self.lastblock:
            raise StopIteration

        if block == self.firstblock:
            start_offset = self.firstblock_start_offset
        else:
            start_offset = 0

        if block == self.lastblock:
            stop_offset = self.lastblock_stop_offset
        else:
            stop_offset = self.blocksize
        # print(f'block: {block}, starto: {start_offset}, stopo {stop_offset}')

        self.i += 1
        return fa[block, start_offset:stop_offset,:][:, self.channels]


class IterData(IterRawData):
    """
    Iterate over blocks of data, scaled with sensitivity and integer scaling
    between 0 and 1
    """

    def __init__(self, fa, channels, sensitivity, **kwargs):
        super().__init__(fa, channels, **kwargs)
        self.sens = np.asarray(sensitivity)[self.channels]
        assert self.sens.ndim == 1

    def __next__(self):
        nextraw = super().__next__()
        return scaleBlockSens(nextraw, self.sens)


class Measurement:
    """Provides access to measurement data stored in the h5 measurement file
    format."""

    # Store a dict of open measurements, with uuid string as a key. We store them as a weak ref.
    uuid_s = WeakValueDictionary()

    def __init__(self, fn):
        """Initialize a Measurement object based on the filename."""

        # Add extension if tried to open without exension
        if DOTMEXT not in fn:
            fn += DOTMEXT

        # Full filepath
        self.fn = fn

        # Open the h5 file in read-plus mode, to allow for changing the
        # measurement comment.
        with h5.File(fn, "r") as f:
            # Check for video data
            try:
                f["video"]
                self.has_video = True
            except KeyError:
                self.has_video = False

            self.nblocks, self.blocksize, self.nchannels = f["audio"].shape
            dtype = f["audio"].dtype
            self.dtype = dtype
            self.sampwidth = getSampWidth(dtype)

            self.samplerate = f.attrs["samplerate"]
            self.N = self.nblocks * self.blocksize
            self.T = self.N / self.samplerate

            try:
                self.version_major = f.attrs["LASP_VERSION_MAJOR"]
                self.version_minor = f.attrs["LASP_VERSION_MINOR"]
            except KeyError:
                # No version information stored 
                self.version_major = 0
                self.version_minor = 1

            try:
                # Try to catch UUID (Universally Unique IDentifier)
                self._UUID = f.attrs['UUID']
                # Flag indicating we have to add a new UUID
                create_new_uuid = False
            except KeyError:
                create_new_uuid = True

            try:
                # UUID of the reference measurement. Should be stored as 
                # a lists of tuples, where each tuple is a combination of (<MeasurementType.value>, <uuid_string>, <last_filename>).

                # The last filename is a filename that *probably* is the reference measurement with
                # given UUID. If it is not, we will search for it in the same directory as `this` measurement.
                # If we cannot find it there, we will give up, and remove the field corresponding to this reference measurement type.
                refMeas_list = f.attrs['refMeas']

                # Build a tuple string from it
                self._refMeas = {}
                for (key, val, name) in refMeas_list:
                    self._refMeas[MeasurementType(int(key))] = (val, name)
            except KeyError:
                self._refMeas = {}

            try:
                self._type_int = f.attrs['type_int']
            except KeyError:
                self._type_int = 0

            # Due to a previous bug, the channel names were not stored
            # consistently, i.e. as 'channel_names' and later camelcase.
            try:
                self._channelNames = f.attrs["channelNames"]
            except KeyError:
                try:
                    self._channelNames = f.attrs["channel_names"]
                    logging.info(
                        "Measurement file obtained which stores channel names with *old* attribute 'channel_names'"
                    )
                except KeyError:
                    # No channel names found in measurement file
                    logging.info("No channel name data found in measurement")
                    self._channelNames = [f"Unnamed {i}" for i in range(self.nchannels)]

            # comment = read-write thing
            if "comment" in f.attrs:
                self._comment = f.attrs["comment"]
            else:
                self._comment = ""

            # Sensitivity
            try:
                sens = f.attrs["sensitivity"]
                self._sens = (
                    sens * np.ones(self.nchannels) if isinstance(sens, float) else sens
                )
            except KeyError:
                self._sens = np.ones(self.nchannels)

            # The time is cached AND ALWAYS ASSUMED TO BE AN IMMUTABLE OBJECT.
            # It is also cached. Changing the measurement timestamp should not
            # be done.
            self._time = f.attrs["time"]

            # Quantity stored as channel.
            self._qtys = None

            try:
                qtys_enum_idx = f.attrs["qtys_enum_idx"]
                self._qtys = [SIQtys.fromInt(idx) for idx in qtys_enum_idx]
            except KeyError:
                try:
                    qtys_json = f.attrs["qtys"]
                    # Load quantity data
                    self._qtys = [Qty.from_json(qty_json) for qty_json in qtys_json]
                except KeyError:
                    # If quantity data is not available, this is an 'old'
                    # measurement file.
                    pass

            if self._qtys is None:
                self._qtys = [SIQtys.default() for i in range(self.nchannels)]
                logging.debug(
                    f"Physical quantity data not available in measurement file. Assuming {SIQtys.default}"
                )

        if create_new_uuid:
            # Create and store a random UUID based on *now* and store it forever
            # inside of the Measurement file
            self.genNewUUID()
        else:
            if self.UUID in Measurement.uuid_s.keys():
                raise RuntimeError(f"Measurement '{self.name}' is already opened. Cannot open measurement twice. Note: this error can happen when measurements are manually copied.")

        # Store weak reference to 'self' in list of UUID's. They are removed when no file is open anymore
        Measurement.uuid_s[self._UUID] = self

    def rename(self, newname: str):
        """
        Try to rename the measurement file.

        Args:
            newname: New name, with or without extension
        """
        _ , ext = os.path.splitext(newname)
        # Add proper extension if new name is given without extension.
        if ext != DOTMEXT:
            newname = newname + DOTMEXT
            
        # Folder, Base filename + extension
        folder, _ = os.path.split(self.fn)

        newname_full = str(pathlib.Path(folder) / newname)
        os.rename(self.fn, newname_full)

        self.fn = newname_full

    def genNewUUID(self):
        """
        Create new UUID for measurement and store in file.
        """
        self.setAttribute('UUID', str(uuid.uuid1()))

    @property
    def UUID(self):
        """
        Universally unique identifier
        """
        return self._UUID
    
    def getRefMeas(self, mtype: MeasurementType):
        """
        Return corresponding reference measurement, if configured and can be found. If the reference 
        measurement is currently not open, it tries to open it by traversing other measurement
        files in the current directory. Throws a runtime error in case the reference measurement cannot be found.
        Throws a ValueError when the reference measurement is not configured.
        """
        
        # See if we can find the UUID for the required measurement type
        try:
            required_uuid, possible_name = self._refMeas[mtype]
        except KeyError:
            raise ValueError(f"No reference measurement configured for '{self.name}'")
            
        m = None
        # Try to find it in the dictionary of of open measurements
        if required_uuid in Measurement.uuid_s.keys():
            m = Measurement.uuid_s[required_uuid]
            logging.debug(f'Returned reference measurement {m.name} from list of open measurements')
        
        # Not found in list of openend measurements. See if we can open it using its last stored file name we know of
        if m is None:
            try:
                m = Measurement(possible_name)
                if m.UUID == required_uuid:
                    logging.info(f'Opened reference measurement {m.name} by name')
            except Exception as e:
                logging.error(f'Could not find reference measurement using file name: {possible_name}')

        # Last resort, see if we can find the right measurement in the same folder
        if m is None:
            try:
                folder, _ = os.path.split(self.fn)
                m = Measurement.fromFolderWithUUID(required_uuid, folder, skip=[self.name])
                logging.info('Found reference measurement in folder with correct UUID. Updating name of reference measurement')
                # Update the measurement file name in the list, such that next time it
                # can be opened just by its name.
                self.setRefMeas(m)
            except:
                logging.error("Could not find the reference measurement. Is it deleted?")

        # Well, we found it. Now make sure the reference measurement actually has the right type (User could have marked it as a NotSpecific for example in the mean time).
        if m is not None:
            if m.measurementType() != mtype:
                m.removeRefMeas(mtype)
                raise RuntimeError(f"Reference measurement for {self.name} is not a proper reference (anymore).")
            
            # Whow, we passed all security checks, here we go!
            return m
        else:
            # Nope, not there.
            raise RuntimeError(f"Could not find the reference measurement for '{self.name}'. Is it deleted?")

    def removeRefMeas(self, mtype: MeasurementType):
        """
        Remove an existing reference measurement of specified type from this measurement. Silently ignores this
        action if no reference measurement of this type is configured.

        """
        try:
            del self._refMeas[mtype]
            self.__storeReafMeas()
        except KeyError:
            pass

    def __storeReafMeas(self):
        """
        Internal method that syncs the dictionary of reference methods with the backing HDF5 file
        """
        with self.file("r+") as f:
            # Update attribute in file. Stored as a flat list of string tuples:
            # [(ref_value1, uuid_1, name_1), (ref_value2, uuid_2, name_2), ...]
            reflist = list((str(key.value), val1, val2) for key, (val1, val2) in self._refMeas.items())
            # print(reflist)
            f.attrs['refMeas'] = reflist

    def setRefMeas(self, m: Measurement):
        """
        Set a reference measurement for the given measurement. If this measurement is already
        a reference measurement, the previous reference measurement type is overwritten, such that there is
        only one measurement that is the reference of a certain 'MeasurementType'
        """
        mtype = m.measurementType()
        if mtype == MeasurementType.NotSpecific:
            raise ValueError('Measurement to be set as reference is not a reference measurement')
    
        self._refMeas[mtype] = (m.UUID, m.name)
        self.__storeReafMeas()

    @staticmethod
    def fromFolderWithUUID(uuid_str: str, folder: str='', skip=[]):
        """
        Returns Measurement object from a given UUID string. It first tries to find whether there
        is an uuid in the static list of weak references. If not, it will try to open files in 
        the current file path.
        """
        for fn in glob.glob(str(pathlib.Path(folder)) + f'/*{DOTMEXT}'):
            # Do not try to open this file in case it is in the 'skip' list.
            if len(list(filter(lambda a: a in fn, skip))) > 0:
                continue
                
            try:
                m = Measurement(fn)
                if m.UUID == uuid_str: 
                    # Update 'last_fn' attribute in dict of stored reference measurements
                    return m
            except Exception as e:
                logging.error(f'Possible measurement file {fn} returned error {e} when opening.')

        raise RuntimeError(f'Measurement with UUID {uuid_str} could not be found.')

    def setAttribute(self, attrname: str, value):
        """
        Set an attribute in the measurement file, and keep a local copy in
        memory for efficient accessing.

        Args:
            attrname: name of attribute, a string
            value: the value. Should be anything that can be stored as an attribute in HDF5.
        """
        with self.file("r+") as f:
            # Update comment attribute in the file
            f.attrs[attrname] = value
        setattr(self, "_" + attrname, value)

    def isType(self, type_: MeasurementType):
        """
        Returns True when a measurement is flagged as being of a certaint "MeasurementType"
        """
        if (type_.value & self._type_int):
            return True
        elif type_.value == self._type_int == 0:
            return True
        return False
    
    def setType(self, type_: MeasurementType):
        """
        Set the measurement type to given type
        """
        self.setAttribute('type_int', type_.value)

    def measurementType(self):
        """
        Returns type of measurement
        """
        return MeasurementType(self._type_int)

    @property
    def name(self):
        """Returns filename base without extension."""
        _, fn = os.path.split(self.fn)
        return os.path.splitext(fn)[0]

    @property
    def channelNames(self):
        return self._channelNames

    @channelNames.setter
    def channelNames(self, newchnames):
        """
        Returns list of the names of the channels
        """
        if len(newchnames) != self.nchannels:
            raise RuntimeError("Invalid length of new channel names")
        self.setAttribute("channelNames", newchnames)

    @property
    def channelConfig(self):
        """
        Returns list of current channel configuration data.
        """
        chcfg = []
        for chname, sens, qty in zip(self.channelNames, self.sensitivity, self.qtys):
            ch = DaqChannel()
            ch.enabled = True
            ch.name = chname
            ch.sensitivity = sens
            ch.qty = qty.cpp_enum
            chcfg.append(ch)
        return chcfg

    @channelConfig.setter
    def channelConfig(self, chcfg: List[DaqChannel]):
        """
        Set new channel configuration from list of DaqChannel objects. 

        Use cases:
            - Update channel types, sensitivities etc.

        Args:
            chchfg: New channel configuration
        """
        if len(chcfg) != self.nchannels:
            raise RuntimeError("Invalid number of channels")

        chname = []
        sens = []
        qtys = []
        for ch in chcfg:
            chname.append(ch.name)
            sens.append(ch.sensitivity)
            qtys.append(SIQtys.fromCppEnum(ch.qty))

        self.channelNames = chname
        self.sensitivity = sens
        self.qtys = qtys

    @property
    def qtys(self):
        return self._qtys

    @qtys.setter
    def qtys(self, newqtys):
        if not len(newqtys) == len(self._qtys):
            raise ValueError("Invalid number of quantities")
        qtys_int = [qty.toInt() for qty in newqtys]
        # Use setAttribute here, but thos store the jsonified version as well,
        # which we have to overwrite again with the deserialized ones. This is
        # actually not a very nice way of coding.
        with self.file("r+") as f:
            # Update comment attribute in the file
            f.attrs["qtys_enum_idx"] = qtys_int

        self._qtys = newqtys

    @contextmanager
    def file(self, mode="r"):
        """Contextmanager which opens the storage file and yields the file.

        Args:
            mode: Opening mode for the file. Should either be 'r', or 'r+'
        """
        if mode not in ("r", "r+"):
            raise ValueError("Invalid file opening mode.")
        with h5.File(self.fn, mode) as f:
            yield f

    @property
    def comment(self):
        """Return the measurement comment.

        Returns:
            The measurement comment (text string)
        """
        return self._comment

    @comment.setter
    def comment(self, cmt):
        """Set the measurement comment.

        Args:
            cmt: Comment text string to set
        """
        with self.file("r+") as f:
            # Update comment attribute in the file
            f.attrs["comment"] = cmt
            self._comment = cmt

    @property
    @lru_cache()
    def recTime(self):
        """Returns the total recording time of the measurement, in float
        seconds."""
        return self.blocksize * self.nblocks / self.samplerate

    @property
    def time(self):
        """Returns the measurement time in seconds since the epoch."""
        return self._time

    @property
    @lru_cache()
    def timestr(self):
        """
        Return a properly formatted string of the measurement time, in order of

        year-month-day hour etc.

        """
        time_struct = time.localtime(self.time)
        time_string = time.strftime("%Y-%m-%d %H:%M:%S", time_struct)
        return time_string

    def rms(self, channels=None, substract_average=False):
        """Returns the root mean square values for each channel

        Args:
            channels: list of channels
            substract_average: If set to true, computes the rms of only the
            oscillating component of the signal, which is in fact the signal
            variance.

        Returns:
            1D array with rms values for each channel
        """
        meansquare = 0.0  # Mean square of the signal, including its average
        sum_ = 0.0  # Sumf of the values of the signal, used to compute average
        N = 0
        with self.file() as f:
            for block in self.iterData(channels):
                Nblock = block.shape[0]
                sum_ += np.sum(block, axis=0)
                N += Nblock
                meansquare += np.sum(block ** 2, axis=0) / self.N

        avg = sum_ / N
        # In fact, this is not the complete RMS, as in includes the DC
        # If p = p_dc + p_osc, then rms(p_osc) = sqrt(ms(p)-ms(p_dc))
        if substract_average:
            meansquare -= avg ** 2
        rms = np.sqrt(meansquare)
        return rms

    def variance(self, channels=None):
        return self.rms(substract_average=True)

    def rawData(self, channels=None, **kwargs):
        """Returns the raw data as stored in the measurement file,
        without any transformations applied

        args:
            channels: list, or tuple of channel numbers to export. If not defined, all
            channels in the measurement are returned

        returns:
            Numpy array with data. The first axis is always the time instance,
            the second axis the channel number.

        """
        if channels is None:
            channels = list(range(self.nchannels))

        rawdata = []

        with self.file() as f:
            for block in IterRawData(f, channels, **kwargs):
                rawdata.append(block)

        return np.concatenate(rawdata, axis=0)

    def iterData(self, channels, **kwargs):
        sensitivity = kwargs.pop("sensitivity", self.sensitivity)
        if channels is None:
            channels = list(range(self.nchannels))
        with self.file() as f:
            for block in IterData(f, channels, sensitivity, **kwargs):
                yield block

    def data(self, channels=None, **kwargs):
        """
        Returns the measurement data, scaled and sensitivity applied.
        """
        data = []
        for d in self.iterData(channels, **kwargs):
            data.append(d)

        return np.concatenate(data, axis=0)

    def CPS(self, channels=None, **kwargs):
        """
        Compute single-sided Cross-Power-Spectrum of the measurement channels

        Args:
            channels: Channels to compute for (numbers)

        Optional arguments:
            nfft: FFT length
            window: Window type
            overlap: Overlap percentage (value between 0.0 and up to and
            including 100.0)
            weighting:

        Returns:
            Cross-power-spectra. C[freq, ch_i, ch_j] = C_ij

        """
        nfft = kwargs.pop("nfft", 2048)
        window = kwargs.pop("windowType", Window.WindowType.Hann)
        overlap = kwargs.pop("overlap", 50.0)

        if channels is None:
            channels = list(range(self.nchannels))

        nchannels = len(channels)
        aps = AvPowerSpectra(nfft, window, overlap)
        freq = getFreq(self.samplerate, nfft)

        for data in self.iterData(channels, **kwargs):
            CS = aps.compute(data)

        return freq, aps.get_est()

    def periodicAverage(self, N, channels=None, noiseCorrection=True, **kwargs):
        """
        Return the (coherent) periodic average the measurement. This method is
        useful for a situation of periodic excitation.

        Args:
            N: The number of samples in one period. This value should
            correspond with the period of the excitation!
            noiseCorrection: whether to apply coherent averaging, according to
            the Sliding Window correlation method (SWiC): Telle et al.: A Novel
            Approach for Impulse Response Measurements with Time-Varying Noise.
            If set to False, just the arithmetic average is used.
        """
        # Create blocks of equal length N
        Ntot = self.N
        Nblocks = Ntot // N

        # TODO: This method graps the whole measurement file into memory. Can
        # only be done with relatively small measurement files.
        signal = self.data(channels)

        # Estimate noise power in block
        blocks = [signal[i * N: (i + 1) * N] for i in range(Nblocks)]

        if noiseCorrection:
            # The difference between the measured signal in the previous block and
            # the current block
            en = [None] + [blocks[i] - blocks[i - 1] for i in range(1, Nblocks)]

            noise_est = [None] + [
                -np.average(en[i] * en[i + 1]) for i in range(1, len(en) - 1)
            ]

            # Create weighting coefficients
            sum_inverse_noise = sum([1 / n for n in noise_est[1:]])
            c_n = [1 / (ni * sum_inverse_noise) for ni in noise_est[1:]]
        else:
            c_n = [1 / (Nblocks - 2)] * (Nblocks - 2)

        assert np.isclose(sum(c_n), 1.0)
        assert Nblocks - 2 == len(c_n)

        # Average signal over blocks
        avg = np.zeros((blocks[0].shape), dtype=float)
        for n in range(0, Nblocks - 2):
            avg += c_n[n] * blocks[n + 1]

        return avg

    def periodicCPS(self, N, channels=None, **kwargs):
        """
        Compute Cross-Spectral Density based on periodic excitation. Uses noise
        reduction by time-averaging the data.
        """

        if channels is None:
            channels = list(range(self.nchannels))

        nchannels = len(channels)
        window = Window.rectangular
        ps = PowerSpectra(N, window)

        avg = np.asfortranarray(self.periodicAverage(N, channels, **kwargs))
        CS = ps.compute(avg)
        freq = getFreq(self.samplerate, N)

        return freq, CS

    @property
    def sensitivity(self):
        """Sensitivity of the data in U^-1, from floating point data scaled
        between -1.0 and 1.0 to Units [U].

        If the sensitivity is not stored in the measurement file, this
        function returns 1.0 for each channel
        """
        return self._sens

    @sensitivity.setter
    def sensitivity(self, sens):
        """Set the sensitivity of the measurement in the file.

        Args:
            sens: sensitivity data, should be a float, or an array of floats
                  equal to the number of channels.
        """
        if isinstance(sens, float):
            # Put all sensitivities equal
            sens = sens * np.ones(self.nchannels)
        elif isinstance(sens, list):
            sens = np.asarray(sens)

        valid = sens.ndim == 1
        valid &= sens.shape[0] == self.nchannels
        valid &= sens.dtype == float
        if not valid:
            raise ValueError("Invalid sensitivity value(s) given")
        with self.file("r+") as f:
            f.attrs["sensitivity"] = sens
        self._sens = sens

    def checkOverflow(self, channels):
        """Coarse check for overflow in measurement.

        Return:
            True if overflow is possible, else False
        """

        for block in self.iterData(channels):
            dtype = block.dtype
            if dtype.kind == "i":
                # minvalue = np.iinfo(dtype).min
                maxvalue = np.iinfo(dtype).max
                if np.max(np.abs(block)) >= 0.9 * maxvalue:
                    return True
            else:
                # Cannot check for floating point values.
                return False
            return False

    def exportAsWave(self, fn=None, force=False, dtype=None, normalize=False, **kwargs):
        """Export measurement file as wave. In case the measurement data is
        stored as floats, the values are scaled to the proper integer (PCM)
        data format.

        Args:
            fn: If given, this will be the filename to write to. If the
            filename does not end with '.wav', this extension is added.

            force: If True, overwrites any existing files with the given name
            , otherwise a RuntimeError is raised.

            dtype: if not None, convert data to this data type.
            Options are 'int16', 'int32', 'float32'.

            normalize: If set: normalize the level to something sensible.
        """
        if fn is None:
            fn = self.fn
            fn = os.path.splitext(fn)[0]

        if os.path.splitext(fn)[1] != ".wav":
            fn += ".wav"

        if os.path.exists(fn) and not force:
            raise RuntimeError(f"File already exists: {fn}")

        if not np.isclose(self.samplerate % 1, 0):
            raise RuntimeError(
                f"Sample rates should be approximately integer for exporting to Wave to work"
            )

        # TODO: With VERY large measurment files, this is not possible! Is this
        # a theoretical case?
        # TODO: add sensitivity? Then use self.data() instead of self.rawData()
        data = self.rawData(**kwargs)

        if normalize:
            # Scale back to maximum of absolute value
            maxabs = np.max(np.abs(data))
            data = data / maxabs  # "data /= maxabs" fails if dtpyes differ

        if dtype == None:
            dtype = data.dtype  # keep existing
            logging.debug(f"dtype not passed as arg; using dtype = {dtype}")

        # dtype conversion
        if dtype == "int16":
            newtype = np.int16
            newsampwidth = 2
        elif dtype == "int32":
            newtype = np.int32
            newsampwidth = 4
        elif dtype == "float32":
            newtype = np.float32
        elif dtype == "float64":
            newtype = np.float64
        else:
            logging.debug(f"cannot handle this dtype {dtype}")
            pass

        # Convert range to [-1, 1]
        # TODO: this is wrong for float data where full scale > 1
        sensone = np.ones_like(self.sensitivity)
        data = scaleBlockSens(data, sensone)

        if dtype == "int16" or dtype == "int32":
            # Scale data to integer range and convert to integers
            scalefac = 2 ** (8 * newsampwidth - 1) - 1
            data = (data * scalefac).astype(newtype)

        wavfile.write(fn, int(self.samplerate), data.astype(newtype))

    @staticmethod
    def fromFile(fn):
        """
        Try to open measurement from a given file name. First checks
        whether the measurement is already open. Otherwise it might
        happen that a Measurement object is created twice for the same backing file, which we do not allow.
        """
        # See if the base part of the filename is referring to a file that is already open
        with h5.File(fn, 'r') as f:
            try:
                theuuid = f.attrs['UUID']
            except KeyError:
                # No UUID stored in measurement. This is an old measurement that did not have UUID's
                # We create a new UUID here such that the file is opened from the filesystem
                # anyhow.
                theuuid = str(uuid.uuid1())

        if theuuid in Measurement.uuid_s.keys():
            return Measurement.uuid_s[theuuid]
        
        return Measurement(fn)

    @staticmethod
    def fromtxt(
        fn,
        skiprows,
        samplerate,
        sensitivity,
        mfn=None,
        timestamp=None,
        delimiter="\t",
        firstcoltime=True,
    ):
        """Converts a txt file to a LASP Measurement file, opens the associated
        Measurement object and returns it. The measurement file will have the
        same file name as the txt file, except with h5 extension.

        Args:
            fn: Filename of text file
            skiprows: Number of header rows in text file to skip
            samplerate: Sampling frequency in [Hz]
            sensitivity: 1D array of channel sensitivities
            mfn: Filepath where measurement file is stored. If not given,
            a h5 file will be created along fn, which shares its basename
            timestamp: If given, a custom timestamp for the measurement
            (integer containing seconds since epoch). If not given, the
            timestamp is obtained from the last modification time.
            delimiter: Column delimiter
            firstcoltime: If true, the first column is the treated as the
            sample time.
        """
        if not os.path.exists(fn):
            raise ValueError(f"File {fn} does not exist.")
        if timestamp is None:
            timestamp = os.path.getmtime(fn)

        if mfn is None:
            mfn = os.path.splitext(fn)[0] + DOTMEXT
        else:
            mfn = os.path.splitext(mfn)[0] + DOTMEXT

        dat = np.loadtxt(fn, skiprows=skiprows, delimiter=delimiter)
        if firstcoltime:
            time = dat[:, 0]
            if not np.isclose(time[1] - time[0], 1 / samplerate):
                raise ValueError(
                    "Samplerate given does not agree with " "samplerate in file"
                )

            # Chop off first column
            dat = dat[:, 1:]
        nchannels = dat.shape[1]
        if nchannels != sensitivity.shape[0]:
            raise ValueError(
                f"Invalid sensitivity length given. Should be: {nchannels}"
            )

        with h5.File(mfn, "w") as hf:
            hf.attrs["samplerate"] = samplerate
            hf.attrs["sensitivity"] = sensitivity
            hf.attrs["time"] = timestamp
            hf.attrs["blocksize"] = 1
            hf.attrs["nchannels"] = nchannels
            ad = hf.create_dataset(
                "audio",
                (1, dat.shape[0], dat.shape[1]),
                dtype=dat.dtype,
                maxshape=(1, dat.shape[0], dat.shape[1]),
                compression="gzip",
            )
            ad[0] = dat
        return Measurement(mfn)

    @staticmethod
    def fromnpy(
        data,
        samplerate,
        sensitivity,
        mfn,
        timestamp=None,
        qtys: List[SIQtys]=None,
        channelNames: List[str]=None,
        force=False,
    ) -> Measurement:
        """
        Converts a numpy array to a LASP Measurement file, opens the
        associated Measurement object and returns it. The measurement file will
        have the same file name as the txt file, except with h5 extension.

        Args:
            data: Numpy array, first column is sample, second is channel. Can
            also be specified with a single column for single-channel data.

            samplerate: Sampling frequency in [Hz]

            sensitivity: 1D array of channel sensitivities in [U^-1], where U is
            the recorded unit.

            mfn: Filepath of the file where the data is stored.

            timestamp: If given, a custom timestamp for the measurement
            (integer containing seconds since epoch).

            qtys: If a list of physical quantity data is given here

            channelNames: Name of the channels

            force: If True, overwrites existing files with specified `mfn`
            name.
        """

        if os.path.splitext(mfn)[1] != DOTMEXT:
            mfn += DOTMEXT
        if os.path.exists(mfn) and not force:
            raise ValueError(f"File {mfn} already exist.")
        if timestamp is None:
            timestamp = int(time.time())

        if data.ndim != 2:
            data = data[:, np.newaxis]

        try:
            len(sensitivity)
        except:
            raise ValueError("Sensitivity should be given as array-like data type")
        sensitivity = np.asarray(sensitivity)

        nchannels = data.shape[1]
        if nchannels != sensitivity.shape[0]:
            raise ValueError(
                f"Invalid sensitivity length given. Should be: {nchannels}"
            )

        if channelNames is not None:
            if len(channelNames) != nchannels:
                raise RuntimeError("Illegal length of channelNames list given")

        if qtys is None:
            qtys = [SIQtys.AP] * nchannels
        else:
            if len(qtys) != nchannels:
                raise RuntimeError("Illegal length of qtys list given")

        with h5.File(mfn, "w") as hf:
            hf.attrs["samplerate"] = samplerate
            hf.attrs["sensitivity"] = sensitivity
            hf.attrs["time"] = timestamp
            hf.attrs["blocksize"] = 1
            hf.attrs["nchannels"] = nchannels

            # Add physical quantity indices
            hf.attrs['qtys_enum_idx'] = [qty.toInt() for qty in qtys]

            # Add channel names in case given
            if channelNames is not None:
                hf.attrs["channelNames"] = channelNames

            ad = hf.create_dataset(
                "audio",
                (1, data.shape[0], data.shape[1]),
                dtype=data.dtype,
                maxshape=(1, data.shape[0], data.shape[1]),
                compression="gzip",
            )
            ad[0] = data
        return Measurement(mfn)

    @staticmethod
    def fromWaveFile(fn, newfn=None, force=False, timestamp=None):
        """Convert a measurement file to a wave file, and return the
        measurement handle."""
        if timestamp is None:
            timestamp = int(time.time())

        base_fn = os.path.splitext(fn)[0]
        if newfn is None:
            newfn = base_fn + ".h5"
        if os.path.exists(newfn) and not force:
            raise RuntimeError(
                f'Measurement file name {newfn} already exists in path, set "force" to true to overwrite'
            )

        samplerate, data = wavfile.read(fn)
        if data.ndim == 2:
            nframes, nchannels = data.shape
        else:
            nchannels = 1
            nframes = len(data)
            data = data[:, np.newaxis]
        sensitivity = np.ones(nchannels)

        with h5.File(newfn, "w") as hf:
            hf.attrs["samplerate"] = samplerate
            hf.attrs["nchannels"] = nchannels
            hf.attrs["time"] = timestamp
            hf.attrs["blocksize"] = 1
            hf.attrs["sensitivity"] = sensitivity
            ad = hf.create_dataset(
                "audio",
                (1, nframes, nchannels),
                dtype=data.dtype,
                maxshape=(1, nframes, nchannels),
                compression="gzip",
            )
            ad[0] = data

        return Measurement(newfn)