Pre-Processing II: Aligning 16 kHz audio with 48 khZ audio from a second source

Overview

The LSL setup records audio at 16 kHz, which is sufficient for speech intelligibility, but not for high-quality audio analysis. Therefore, a second audio recording was made using Audacity, with sampling rate 48 kHz. This notebook aligns the two audio streams based on cross-correlation, allowing us to replace the lower-quality audio from the XDF files with the higher-quality audio from Audacity.

Additionally, we denoise the audio to remove background noise.

Lastly, we merge the audio and video files we create in the previous script to allow for checking of audio-video synchrony and other potential issues that need to be sorted out before further analysis.

Aligning audios with shign

We will now align the 48kHz audio to the 16kHz which has been recorded via LSL and therefore is synced with all remaining streams. We will use shign package (see Github)

To verify the alignment, we can plot parts of the two signals after the procedure. If the alignment is successful, the waveforms should closely match each other. Files that fail the check will be aligned manually in Adobe Premiere Pro.

from __future__ import print_function
import numpy as np
import wave
import soundfile as sf
import noisereduce as nr
from scipy.io import wavfile
import glob
import os
import shign
import matplotlib.pyplot as plt
import wave
import pandas as pd
import gc
from scipy.signal import find_peaks, peak_widths
import struct
import ffmpeg
import noisereduce as nr

curfolder = os.getcwd()
folder16 = os.path.join(curfolder, 'data', 'Data_processed', 'CsvDataTS_raw', 'Audio') 
folder48 = os.path.join(curfolder, '..', "00_raw")

# Collect all 16k files
files16 = glob.glob(os.path.join(folder16, '*.wav'))
files16 = [f for f in files16 if '_denoised' in os.path.basename(f)] # keep only denoised files

# Collect all 48k files
files48 = glob.glob(os.path.join(folder48, '*', '*.wav'))
files48 = [f for f in files48 if '_denoised' not in os.path.basename(f) and '_aligned' not in os.path.basename(f)] # ignore all denoised if already available

print(f'Found {len(files16)} denoised 16k files and {len(files48)} 48k files.')

errorlist = []
skip = ['4_2', '23_2', '15_1', '16_2']
PLOT_DEBUG = True # set this to False if you don't need plots for every file

for file48 in files48:
    print('working on ' + file48)

    # get session ID
    basename = os.path.basename(file48).split('.')[0]
    if '_reinitiated_' in basename:
        sessionID = '_'.join(basename.split('_')[0:3])
        corresponding16 = [f for f in files16 if os.path.basename(f).startswith(sessionID)]
    else:
        sessionID = '_'.join(basename.split('_')[0:2])
        corresponding16 = [f for f in files16
                           if os.path.basename(f).startswith(sessionID)]
        corresponding16 = [f for f in corresponding16 if '_rein_' not in os.path.basename(f)]
        
    if sessionID in skip:
        print('Session ' + sessionID + ' is in skip list, skipping.')
        continue

    if len(corresponding16) == 0:
        print('No corresponding 16k file found for ' + file48)
        errorlist.append('No 16k for ' + file48)
        continue

    try:
        # Find the corresponding 16kHz file
        correspondingfile = corresponding16[0]
        print('corresponding 16k file: ' + correspondingfile)

        # Open wav file
        with wave.open(file48, 'rb') as audio48, wave.open(correspondingfile, 'rb') as audio16:
            print('opening files')

            # Check number of channels
            n_channels48 = audio48.getnchannels()
            n_channels16 = audio16.getnchannels()

            # Sample rates
            sample_rate48 = audio48.getframerate()
            sample_rate16 = audio16.getframerate()

            # Read raw bytes
            raw48 = audio48.readframes(-1)
            raw16 = audio16.readframes(-1)

        # Convert to float32 and normalize to [-1, 1] to avoid overflow
        audio_data48 = np.frombuffer(raw48, dtype=np.int16)
        if n_channels48 > 1:
            audio_data48 = audio_data48.reshape(-1, n_channels48)[:, 0]  # take only first channel
        audio_data48 = audio_data48.astype(np.float32) / 32768.0

        audio_data16 = np.frombuffer(raw16, dtype=np.int16)
        if n_channels16 > 1:
            audio_data16 = audio_data16.reshape(-1, n_channels16)[:, 0]  # take only first channel
        audio_data16 = audio_data16.astype(np.float32) / 32768.0

        # Free raw byte buffers
        del raw48, raw16

        print('aligning')
        my_aligned_audio1, my_aligned_audio2 = shign.shift_align(
            audio_a=audio_data16,
            audio_b=audio_data48,
            sr_a=sample_rate16,
            sr_b=sample_rate48,
            align_how='pad_and_crop_one_to_match_other'
        )

        # Optional plotting
        if PLOT_DEBUG:
            plt.figure(figsize=(20, 10))
            plt.subplot(2, 1, 1)
            plt.plot(my_aligned_audio1)
            plt.title('aligned audio 1')
            plt.subplot(2, 1, 2)
            plt.plot(my_aligned_audio2)
            plt.title('aligned audio 2')
            plt.tight_layout()
            plt.show()
            plt.close()

        # Save the output
        work_dir = os.path.dirname(file48)
        wav2_path = os.path.join(
            work_dir,
            os.path.basename(file48).replace('.wav', '_aligned.wav')
        )
        sf.write(wav2_path, my_aligned_audio2, 48000)

        # Explicitly free big arrays and trigger GC
        del audio_data48, audio_data16, my_aligned_audio1, my_aligned_audio2
        gc.collect()

    except Exception as e:
        print(f'Error in file: {file48}')
        print(f'  -> {type(e).__name__}: {e}')
        errorlist.append('Error for in syncing for ' + file48)
        gc.collect()
        continue

Note that if the alignment is visibly off, we align the two wav files manually in Adobe Premiere using the function Synchronize...

Getting rid of background noise and other artifacts

Now when all audio 48 kHz files are perfectly aligned with remaining (LSL) streams, we need to denoise them to get rid of background noise and other unwanted artifacts.

Upon manual inspection, we found that most of the noise can be attributed to stationary background noise (e.g., air conditioning). We will denoise those files using the noisereduce package, using different strengths depending on the amount of noise present in each file.

Some of the files have been found to contain electrical interference artifacts (i.e., hiss). Those will be treated separately in the next section, using notch filters.

alignedfolder = os.path.join(curfolder, '..', '00_raw')
alignedfiles = glob.glob(os.path.join(alignedfolder, '*', '*_aligned*.wav'))
alignedfiles = [f for f in alignedfiles if '_denoised' not in os.path.basename(f)] # get rid of denoised if already available

print(f'Found {len(alignedfiles)} aligned files.')

Found 146 aligned files.

# Skip all sessions from 3 to 16 - these will go through dehissing pipeline
toskip = ['3_1', '3_2', '4_1', '4_2', '5_1', '5_2', '6_1', '6_2',
          '7_1', '7_2', '8_1', '8_2', '9_1', '9_2', '10_1', '10_2',
          '11_1', '11_2', '12_1', '12_2', '13_1', '13_2', '14_1', '14_2',
          '15_1', '15_2', '16_1', '16_2']

# Here we specify strength of the filter that each session (upon manual auditory inspection) needs
verylittle = ['38_2', '39_1']

little = ['36_2', '37_1', '37_2', '38_1', '39_1',
        '40_2', '41_1', '41_2', '43_1', '43_2', '44_1', '44_2',
        '45_1', '45_2', '46_1', '46_2', '47_1', '47_2', '48_1', '48_2',
        '49_1', '49_2', '50_1', '50_2', '51_1', '51_2', '52_1', '52_2',
        '53_1', '53_2', '54_2', '55_1', '55_2', '56_1', '56_2',
        '57_1', '57_2', '58_1', '58_2', '59_1', '59_2', '60_2',
        '61_1', '61_2', '62_1', '62_2', '63_1', '63_2', '64_1', '64_2',
        '65_1', '65_2', '66_1', '66_2', '67_1', '69_1', '69_2',
        '71_1', '71_2', '72_1', '72_2']

more = ['27_1', '28_1', '28_2', '29_1', '29_2', '30_1', '30_2',
        '31_1', '31_2', '32_1', '33_1', '33_2', '34_1', '34_2',
        '35_1', '35_2', '36_1']

big = ['1_1', '1_2', '2_1', '2_2','54_1', '67_2', '27_2','32_2', '42_1','42_2', '40_1', # nonstanionary
       '17_1', '17_2', '18_1', '18_2', '19_1', '19_2', '20_1', '20_2', 
         '21_1', '21_2', '22_1', '22_2', '23_1', '23_2', '24_1', '24_2',
         '26_1']

for alignedfile in alignedfiles:
    
    basename = os.path.basename(alignedfile).split('.')[0]
    if basename.startswith(tuple(toskip)):
        print(f'Skipping denoising for {alignedfile}')
        continue

    denoised_path = alignedfile.replace('_aligned.wav', '_aligned_denoised.wav')

    if os.path.exists(denoised_path):
        print(f'Denoised file already exists for {alignedfile}, skipping.')
        continue

    print(f'Denoising {alignedfile}')

    data, rate = sf.read(alignedfile)

    # Create 1 second noise sample from end
    noise_sample = data[-int(1 * rate):]      

    # Perform noise reduction
    if basename.startswith(tuple(verylittle)):
        print('Denoising level 0')
        reduced_noise = nr.reduce_noise(y=data, sr=rate, prop_decrease=0.65, n_std_thresh_stationary=1.8, stationary=True)

    elif basename.startswith(tuple(little)):
        print('Denoising level 1')
        reduced_noise = nr.reduce_noise(y=data, sr=rate, prop_decrease=0.85, n_std_thresh_stationary=1.7, stationary=True)

    elif basename.startswith(tuple(more)):
        print('Denoising level 2')
        reduced_noise = nr.reduce_noise(y=data, sr=rate, prop_decrease=0.87, n_std_thresh_stationary=1.5, stationary=True)
        
    elif basename.startswith(tuple(big)):
        print('Denoising level 3')
        reduced_noise = nr.reduce_noise(y=data, sr=rate, prop_decrease=0.9, n_std_thresh_stationary=1.2, stationary=False)

    else:
        print(f'No specific denoising parameters for {alignedfile}.')
        reduced_noise = nr.reduce_noise(y=data, sr=rate, prop_decrease=0.85, n_std_thresh_stationary=1.7, stationary=True)
        
    # Save denoised audio
    sf.write(denoised_path, reduced_noise, rate)

Stronger denoising: Fast-Fourier Transform (FFT) based noise reduction

Very often, gating algorithms such as noisereduce are not sufficient to remove tonal electrical noise (hissing, humming). This is because such noise is often very strong and concentrated at specific frequencies.

For this reason, we will apply Fast Fourier Transform (FFT) to identify the main frequencies of the noise, and then apply notch filters to attenuate those frequencies. To understand FFT, you can watch Steve Brunton’s excellent video on YouTube. The following code was also inspired by this StackOverflow post.

def compute_median_spectrum(file_list, n_fft=65536, max_files=40):
    """
    Compute median magnitude spectrum across a set of files.
    Use this to find persistent (environment) spectral lines.
    """
    mags = []
    sr_ref = None
    for i, f in enumerate(file_list):
        if i >= max_files:
            break
        try:
            data, sr = sf.read(f)
        except Exception:
            continue
        if data.ndim == 2:
            data = data.mean(axis=1)
        if sr_ref is None:
            sr_ref = sr

        # Trim or pad to reasonable length if file is huge (keeps memory bounded)
        if len(data) > 10 * sr_ref:  # use first 10s
            data = data[:10 * sr_ref]
        FFT = np.fft.rfft(data, n=n_fft)
        mags.append(np.abs(FFT))
    if len(mags) == 0:
        raise RuntimeError("No spectra computed")
    mags = np.vstack(mags)
    median_mag = np.median(mags, axis=0)
    freq = np.fft.rfftfreq(n_fft, d=1.0 / sr_ref)
    return freq, median_mag, sr_ref

def detect_persistent_peaks(freq, median_mag, min_freq=1000, prominence_frac=0.05, min_sep_hz=20):
    """
    Detect peaks in the median spectrum that are likely persistent (environmental) spectral lines.

    Parameters:
    - freq: Frequency bins from FFT
    - median_mag: Median magnitude spectrum across files
    - min_freq: Minimum frequency to consider (default: 1000 Hz)
    - prominence_frac: Fraction of max magnitude to use as prominence threshold
    - min_sep_hz: Minimum separation between peaks in Hz

    Returns:
    - List of frequencies (Hz) of detected persistent peaks
    """

    bins_per_hz = len(freq) / (freq.max() if freq.max()>0 else 1.0)
    distance_bins = max(1, int(min_sep_hz * bins_per_hz))
    prominence = max(1e-12, np.max(median_mag) * prominence_frac)
    peaks, _ = find_peaks(median_mag, prominence=prominence, distance=distance_bins)
    peaks = peaks[freq[peaks] >= min_freq]
    return freq[peaks].tolist()

def attenuate_freq_band(FFT, freq, center_hz, width_hz, attenuation=0.9):
    """
    Attenuate a frequency band in the FFT data using a Hanning window.

    Parameters:
    - FFT: Complex FFT data
    - freq: Frequency bins corresponding to FFT
    - center_hz: Center frequency of band to attenuate
    - width_hz: Width of band to attenuate (total width, not half-width)
    - attenuation: Strength of attenuation (0.0 to 1.0)

    Returns:
    - Modified FFT data with specified band attenuated
    """
    idx_band = np.where((freq >= center_hz - width_hz / 2) & (freq <= center_hz + width_hz / 2))[0]
    if idx_band.size == 0:
        return FFT
    win = np.hanning(len(idx_band))
    multiplier = 1.0 - attenuation * win
    FFT[idx_band] = FFT[idx_band] * multiplier
    return FFT

def erase_freq(
    data,
    rate,
    targets=None,
    plot_before=False,
    plot_after=True,
    min_freq=900,
    min_sep_hz=50,
    prominence_frac=0.05,
    persistent_refs=None,
    match_tol_hz=2.0,
    max_width_hz=10.0,
    profile="balanced"   # choose attenuation profile
):
    """
    Global-FFT tonal attenuation with band-dependent notch parameters.

    If targets provided: attenuate those.
    If targets is None: auto-detect peaks in this recording but ONLY attenuate those that:
      - match a persistent_refs list (if provided), OR
      - are very narrow peaks (width < max_width_hz) indicating tonal electrical lines.

    `profile` controls how strong the attenuation is in each band:
      "very_conservative", "conservative", "balanced", "strong", "very_strong"
    """

    FFT_data = np.fft.rfft(data)
    freq = np.fft.rfftfreq(len(data), d=1.0 / rate)
    magnitude = np.abs(FFT_data)

    # we ignore < 300 Hz here 
    BAND_DEFS = [
        ("low",   300.0, 1200.0),
        ("mid",  1200.0, 3500.0),
        ("high", 3500.0, 8000.0),
        ("ultra",8000.0, np.inf),
    ]

    PROFILES = {
        "very_conservative": {
            "low":   (4.0, 0.25),   # width_hz, attenuation
            "mid":   (6.0, 0.30),
            "high":  (10.0, 0.45),
            "ultra": (15.0, 0.55),
        },
        "conservative": {
            "low":   (5.0, 0.30),
            "mid":   (8.0, 0.40),
            "high":  (15.0, 0.60),
            "ultra": (20.0, 0.70),
        },
        "balanced": {
            "low":   (6.0, 0.35),   
            "mid":   (8.0, 0.45),
            "high":  (15.0, 0.75),
            "ultra": (25.0, 0.90),
        },
        "strong": {
            "low":   (8.0, 0.45),
            "mid":   (12.0, 0.60),
            "high":  (25.0, 0.85),
            "ultra": (35.0, 0.95),
        },
        "very_strong": {
            "low":   (10.0, 0.60),  
            "mid":   (15.0, 0.75),
            "high":  (30.0, 0.95),
            "ultra": (40.0, 0.98),
        },
    }

    if profile not in PROFILES:
        raise ValueError(f"Unknown profile '{profile}'. "
                         f"Choose one of: {list(PROFILES.keys())}")
    profile_params = PROFILES[profile]

    # Helper to map hz → band name
    def band_for_hz(hz):
        if hz < 300.0:
            return None  # ignore sub-300 Hz in this pipeline
        for band_name, lo, hi in BAND_DEFS:
            if lo <= hz < hi:
                return band_name
        return "ultra"  # just in case

    # Choose notch params based on freq & profile
    def choose_notch_params(hz):
        band = band_for_hz(hz)
        if band is None:
            return None, None
        base_width, base_att = profile_params[band]
        # Respect global max_width_hz if given
        if max_width_hz is not None:
            local_width = min(base_width, max_width_hz) if band in ("low", "mid") else base_width
        else:
            local_width = base_width
        local_att = base_att
        return local_width, local_att

    # Peak detection 
    bins_per_hz = len(freq) / (rate / 2.0)
    distance_bins = max(1, int(min_sep_hz * bins_per_hz))
    prominence = max(1e-12, np.max(magnitude) * prominence_frac)
    peaks, _ = find_peaks(magnitude, prominence=prominence, distance=distance_bins)
    peaks = peaks[freq[peaks] >= min_freq]

    if len(peaks) > 0:
        results_width = peak_widths(magnitude, peaks, rel_height=0.5)
        widths_bins = results_width[0]
        widths_hz = widths_bins / bins_per_hz
    else:
        widths_hz = np.array([])

    centers_to_attenuate = []

    if targets:
        centers_to_attenuate = [t for t in targets if t >= min_freq]
    else:
        for i, p in enumerate(peaks):
            hz = freq[p]
            w = widths_hz[i] if i < len(widths_hz) else np.inf
            is_narrow = (w <= max_width_hz)
            matches_persistent = False
            if persistent_refs is not None and len(persistent_refs) > 0:
                if np.any(np.isclose(hz, persistent_refs, atol=match_tol_hz)):
                    matches_persistent = True
            if matches_persistent or is_narrow:
                centers_to_attenuate.append(hz)

    # Plot before attenuation
    if plot_before:
        plt.figure(figsize=(10, 4))
        plt.plot(freq, magnitude, color='gray', label='original')
        if len(peaks):
            plt.scatter(freq[peaks], magnitude[peaks],
                        color='red', zorder=5, label='detected peaks')
        if persistent_refs:
            plt.vlines(persistent_refs, ymin=0, ymax=np.max(magnitude),
                       color='orange', alpha=0.6, label='persistent refs')
        if centers_to_attenuate:
            plt.vlines(centers_to_attenuate, ymin=0, ymax=np.max(magnitude),
                       color='blue', alpha=0.7, label='to attenuate')
        plt.title(f"FFT Magnitude (Before Attenuation, profile='{profile}')")
        plt.xlabel("Frequency (Hz)")
        plt.ylabel("Magnitude")
        plt.xlim(0, min(5000, freq.max()))
        plt.legend()
        plt.grid(True)
        plt.show()

    # Apply attenuation
    for hz in centers_to_attenuate:
        local_width, local_att = choose_notch_params(hz)
        if local_width is None:
            continue
        FFT_data = attenuate_freq_band(FFT_data, freq, hz,
                                       width_hz=local_width,
                                       attenuation=local_att)

    # Plot after attenuation
    if plot_after:
        mag_after = np.abs(FFT_data)
        plt.figure(figsize=(10, 4))
        plt.plot(freq, mag_after, color='magenta', label='after')
        if centers_to_attenuate:
            idx_centers = [int(np.argmin(np.abs(freq - h))) for h in centers_to_attenuate]
            plt.scatter(centers_to_attenuate, mag_after[idx_centers],
                        color='blue', zorder=5, label='attenuated centers')
        plt.title(f"FFT Magnitude After Attenuating Peaks (profile='{profile}')")
        plt.xlabel("Frequency (Hz)")
        plt.ylabel("Magnitude")
        plt.xlim(0, min(5000, freq.max()))
        plt.legend()
        plt.grid(True)
        plt.show()

    new_data = np.fft.irfft(FFT_data)
    return new_data

First, we will compute median spectrum across all noisy files to identify the main frequencies of the noise. We leverage the fact that the noise is consistent across all files, while the speech/vocalization content varies. Therefore, by averaging across all files, we suppress file-specific content (here vocalizations) and highlight stable, persistent lines. Afterwards, we will identify the main (most prominent, well-separated) peaks in the averaged spectrum.

# Keep only those in alignedfiles whose basename starts with skip
alignedfiles = [f for f in alignedfiles if os.path.basename(f).split('.')[0] .startswith(tuple(toskip))]

noisyfiles = []
for f in alignedfiles:
    base = os.path.splitext(os.path.basename(f))[0]  # remove .wav
    session = "_".join(base.split("_")[:2])          # take first two parts (e.g., "24_1")

    if session in toskip:
        noisyfiles.append(f)


print(f'Computing median spectrum from {len(noisyfiles)} noisy files for persistent peak detection.')

Computing median spectrum from 26 noisy files for persistent peak detection.

# Compute median spectrum from example noise files
freq_ref, med_mag, sr_ref = compute_median_spectrum(noisyfiles, n_fft=65536, max_files=30)

# Print the results
print(f"Computed median spectrum from {len(noisyfiles)} files at sample rate {sr_ref} Hz")
print(f"Median spectrum has {len(med_mag)} frequency bins up to {freq_ref[-1]:.1f} Hz")

# Detect persistent peaks in the median spectrum
persistent = detect_persistent_peaks(freq_ref, med_mag, min_freq=900, prominence_frac=0.03)
print(f"Detected {len(persistent)} persistent peaks at frequencies: {persistent}")

Computed median spectrum from 26 files at sample rate 48000 Hz
Median spectrum has 32769 frequency bins up to 24000.0 Hz
Detected 26 persistent peaks at frequencies: [900.146484375, 960.205078125, 999.755859375, 1100.09765625, 1200.439453125, 1250.244140625, 1300.048828125, 1400.390625, 1500.0, 2000.244140625, 3000.0, 3999.755859375, 4999.51171875, 5999.267578125, 6999.0234375, 7998.779296875, 8063.232421875, 8999.267578125, 9335.44921875, 9999.755859375, 10081.787109375, 10999.51171875, 12000.0, 12098.876953125, 18665.771484375, 19995.849609375]

To attenuate the identified noise frequencies, we will apply notch filters at those frequencies. A notch filter is a type of band-stop filter that attenuates a narrow frequency band while leaving other frequencies relatively unaffected. By applying notch filters at the identified noise frequencies, we can effectively reduce the tonal noise while preserving the overall structure of the audio signal.

The erase_freq() function takes in data and compute Fast Fourier Transform (FFT) to convert the signal to frequency domain. It then find peaks in the recording and builds a list of frequencies to attenuate. This list can be either given (targets) or if targets==None, it picks those peaks that are near persistent ones (within match_tol_hz), or narrow enough (max_width_hz). It then applies notch filters at those frequencies by building Hanning windows and scaling the FFT bins accordingly. As a result of Hanning window, the middle of the band is attenuated the most, edges are tapered smoothly. This is to prevent ringing artifacts that would occur with hard brick notch (setting band to zero).

Finally, it converts the signal back to time domain using inverse FFT and also applies a mild noise reduction using noisereduce to clean up any residual noise.

# We will create to profiles, one very conservative and one balanced.
profiles = [
    "very_conservative",
    "balanced",
]

filestoprocess = alignedfiles

for sample in filestoprocess:
    
    print("Processing file:", sample)
    basename = os.path.basename(sample).split('.')[0]

    if not basename.startswith(tuple(toskip)):
        continue
    
    data, rate = sf.read(sample)

    for prof in profiles:
        print("Processing with profile:", prof)
        out = erase_freq(
            data,
            rate,
            min_freq=900,
            min_sep_hz=50,
            prominence_frac=0.05,
            persistent_refs=persistent,   
            match_tol_hz=2.0,
            max_width_hz=10.0,    # caps low/mid band widths
            profile=prof,
            plot_before=True,
            plot_after=True
        )

        # then noisereduce
        reduced = nr.reduce_noise(
            y=out,
            sr=rate,
            stationary=True,
            n_std_thresh_stationary=1.5,
            prop_decrease=0.9,
        )

        out_path = f"{os.path.splitext(sample)[0]}_denoised_{prof}.wav"
        sf.write(out_path, reduced, rate)
        print("Saved", out_path)

Cut 48 kHz session audio to trials

We will now cut the aligned, denoised 48 kHz audio files into trials based on the timestamps extracted from the XDF files in the previous notebook. The resulting trial-wise audio files will be saved for further analysis.

# Set folders
curfolder = os.getcwd()
datafolder = os.path.join(curfolder, 'data')
experiment_to_process = os.path.join(curfolder, '..', '00_raw')

# Collect all 48k denoised files
audio_48_files = glob.glob(os.path.join(experiment_to_process, '**', '*denoised*.wav'), recursive=True)

# This is where we store session 16 kHz audio files
ts_audio = os.path.join(curfolder, 'data', 'Data_processed', 'CsvDataTS_raw')
csvfiles = glob.glob(os.path.join(ts_audio, '*.csv'))
ts_audio_files = [x for x in csvfiles if 'Mic' in x] # keep only those that have Mic

trials = os.path.join(curfolder, 'data', 'Data_processed', 'Data_trials')
csvfiles = glob.glob(os.path.join(trials, '*.csv'))
ts_audio_trials = [x for x in csvfiles if 'Mic' in x] # keep only those that have Mic

# Here we will store the new audio
audio_trials_48 = os.path.join(curfolder, 'data', 'Data_processed', 'Data_trials', 'Audio_48')
if not os.path.exists(audio_trials_48):
    os.makedirs(audio_trials_48)

def find_closest_rows(df, target_number):
    """
    Finds the row(s) in a DataFrame whose first column value is closest to a target number.

    Parameters:
    - df (pd.DataFrame): The DataFrame to search.
    - target_number (float): The number to find the closest value to.

    Returns:
    - pd.DataFrame: A DataFrame containing the row(s) with the closest value to the target number.
    """
    # Calculate the absolute difference between each value in the 'trial_start' column and the target number
    df['abs_difference'] = abs(df.iloc[:, 0] - target_number)
    
    # Find the row with the smallest absolute difference for each row
    closest_rows = df[df['abs_difference'] == df['abs_difference'].min()]

    return closest_rows

# Audio write function
def to_audio(fileloc, timeseries, samplerate = 16000, channels = 1):
    """
    Writes a time series to an audio file.

    Parameters:
    - fileloc (str): The path to the output audio file.
    - timeseries (list or np.array): The time series data to be written to the audio file.
    - samplerate (int, optional): The sample rate of the audio file. Default is 16000.
    - channels (int, optional): The number of channels in the audio file. Default is 1 (mono).
    """
    obj = wave.open(fileloc,'w')
    obj.setnchannels(channels) # mono
    obj.setsampwidth(2)
    obj.setframerate(float(samplerate))
    for i in timeseries:
        data = struct.pack('<h', int(i))
        obj.writeframesraw( data )
    obj.close()

error_list = []
convert = 3 # 48000/16000
no48 = ['15_1', '16_2']

for file in ts_audio_files:
    print('working on ' + file)

    # Get session ID
    if 'rein' in file:
        sessionID = file.split(os.sep)[-1].split('.')[0].split('_')[0] + '_' + file.split(os.sep)[-1].split('.')[0].split('_')[1] + '_rein'
    else:
        sessionID = file.split(os.sep)[-1].split('.')[0].split('_')[0] + '_' + file.split(os.sep)[-1].split('.')[0].split('_')[1]

    if sessionID in no48:
        print('sessionID ' + sessionID + ' has no 48kHz file, skipping...')
        continue

    # Find the corresponding 48kHz file
    try:
        if 'rein' in file:
            files48 = [x for x in audio_48_files if sessionID == (x.split(os.sep)[-1].split('.')[0].split('_')[0] + '_' + x.split(os.sep)[-1].split('.')[0].split('_')[1] + '_rein')]
            # Get rid of those that do NOT have rein in them
            files48 = [x for x in files48 if '_rein_' in x]
            # If there is file that contains 'balanced' take it
            if any('balanced' in x for x in files48):
                file48 = [x for x in files48 if 'balanced' in x][0]
            else:
                file48 = [x for x in files48 if 'denoised' in x][0]

        else:
            file48 = [x for x in audio_48_files if sessionID == (x.split(os.sep)[-1].split('.')[0].split('_')[0] + '_' + x.split(os.sep)[-1].split('.')[0].split('_')[1])]
            file48 = [x for x in file48 if '_rein_' not in x]
            if any('balanced' in x for x in file48):
                file48 = [x for x in file48 if 'balanced' in x][0]
            else:
                file48 = [x for x in file48 if 'denoised' in x][0]

        print('corresponding file: ' + file48)

    except:
        print('no corresponding file found')
        error_list.append('No 48kHz file for ' + sessionID)
        continue

    # Open audio and csv file
    audio48 = wave.open(file48, 'rb')
    data48 = np.frombuffer(audio48.readframes(-1), dtype=np.int16)
    ts_audio = pd.read_csv(file)

    # Get all the audiotrials_ids that have exactly identical sessionID
    if 'rein' in file:
        audiotrials = [
            x for x in ts_audio_trials 
            if sessionID == (x.split(os.sep)[-1].split('.')[0].split('_')[0] + '_' + x.split(os.sep)[-1].split('.')[0].split('_')[1] + '_rein')]
        # Get rid of those that do NOT have rein in them
        audiotrials = [x for x in audiotrials if '_rein_' in x]

    else:
        audiotrials = [
        x for x in ts_audio_trials 
        if sessionID == (x.split(os.sep)[-1].split('.')[0].split('_')[0] + '_' + x.split(os.sep)[-1].split('.')[0].split('_')[1])]
        # Get rid of the one with _rein_
        audiotrials = [x for x in audiotrials if '_rein_' not in x]

    for trial in audiotrials:
        print('working on trial: ' + trial)

        filename_new = trial.split(os.sep)[-1].split('.')[0]
        filename_new = filename_new.replace('16000', '48000')
        file_path = os.path.join(audio_trials_48, filename_new + '.csv')

        wavloc = os.path.join(audio_trials_48, filename_new + '.wav')

        # If it exists, skip
        if os.path.exists(wavloc):
            print(wavloc + ' already exists, skipping...')
            continue

        # Open the csv file
        df = pd.read_csv(trial)

        start = df.iloc[0, 0]
        end = df.iloc[len(df)-1, 0]

        # Find the closest row in ts_audio
        closest_row_start = find_closest_rows(ts_audio, start)
        start_index = closest_row_start.index[0]    
        #print('start index ', start_index)
        closest_row_end = find_closest_rows(ts_audio, end)
        end_index = closest_row_end.index[0]
        #print('end index ', end_index)

        # Convert
        start_index_2 = start_index*convert
        end_index_2 = end_index*convert
        #print('new indices ', start_index_2, end_index_2)

        # Cut the data48
        data2_trial = data48[int(start_index_2):int(end_index_2)]
        
        # Save the wavfile
        to_audio(wavloc, data2_trial, samplerate = 48000)

    # Close the audio and clean the buffer
    audio48.close()
    audio48 = None
    data48 = None
    ts_audio = pd.DataFrame()

# Save error list
if len(error_list) > 0:
    with open(os.path.join(curfolder, 'error_list_cutting_48k.txt'), 'w') as f:
        for item in error_list:
            f.write("%s\n" % item)

Merging audio and video files

Finally, we will merge the denoised, aligned audio files with the corresponding video files created in the previous notebook. This will allow us to review the synchronized audio-visual data for quality assurance and further analysis. We will use ffmpeg for this purpose.

# Set folders
curfolder = os.getcwd()
datafolder = os.path.join(curfolder, 'data')
trialfolder = os.path.join(datafolder, 'Data_processed', 'Data_trials')
videofolder = trialfolder

audiovideo = os.path.join(datafolder, 'Data_processed', 'AudioVideo_48')
if not os.path.exists(audiovideo):
    os.makedirs(audiovideo)

audio_trials_48 = os.path.join(trialfolder, 'Audio_48')
wavfiles48 = glob.glob(os.path.join(audio_trials_48, '*.wav'))

# Inititate error log
errorlog = []

# These sessions will be concatenated with 16kHz because we failed to synchronize the 48 kHz audio with the rest of the streams
with16 = ['15_1', '16_2']

for file in wavfiles48:
    print('Now processing file '+file)

    filename = file.split(os.sep)[-1].split('.')[0]
    dyadIndex = filename.split('_')[0]   # this is dyad number
    partIndex = filename.split('_')[1]

    if 'rein' in filename:
        sessionIndex = dyadIndex + '_' + partIndex + '_rein'
    else:
        sessionIndex = dyadIndex + '_' + partIndex

    if sessionIndex in with16:
        continue

    if '_rein_' in filename: 
        trialIndex = filename.split('_')[4] # this is trial number
        participant = filename.split('_')[8]
        word = filename.split('_')[9]
        modality = filename.split('_')[10]
        
    else:
        trialIndex = filename.split('_')[3] # this is trial number
        participant = filename.split('_')[7] # this is participant 0/1
        word = filename.split('_')[8] # this is the word
        modality = filename.split('_')[9] # this is the modality

    # Handle correction
    if 'c0' in file:
        correction = '_c0'
    elif 'c1' in file:
        correction = '_c1'
    elif 'c2' in file:
        correction = '_c2'
    else:
        correction = ''
    # Trial type
    if '_pr_' in file:
        trialtype = 'pr'
    else:
        trialtype = 'trial'

    # If processing already files that were cut (corrected for faulty start or end)
    if 'cut' in file:
        cut = filename.split('_')[-1]
        cut = '_' + cut
    else:
        cut = ''

    # This is the location where the video will be saved
    output_path = os.path.abspath(os.path.join(audiovideo, f"{sessionIndex}_{trialtype}_{trialIndex}_{participant}_{word}_{modality}{correction}_final.avi"))
  
    if os.path.exists(output_path):
        print('Video already exists, skipping...')
        continue

    # Load in the audio
    print('Loading the audio')
    audio_path = file
    if not os.path.exists(audio_path):
        print(f"Audio file not found: {audio_path}")

    # Input the video with ffmpg
    input_audio = ffmpeg.input(audio_path)

    # Load in the video with matching trialIndex and SessionIndex
    print('Loading the video')
    video_path = os.path.join(videofolder, f"{sessionIndex}_{trialtype}_{trialIndex}_{participant}_{word}_{modality}{correction}{cut}_video_raw.avi")

    if not os.path.exists(video_path):
        print(f"Video file not found: {video_path}")
        errorlog.append('Video file not found: ' + video_path)

    input_video = ffmpeg.input(video_path)

    print('Combining audio and video')

    try:
        ffmpeg.concat(input_video, input_audio, v=1, a=1).output(
            output_path,
            vcodec='libx264',
            acodec='aac',
            video_bitrate='2M',         
            
            ).run(overwrite_output=True)
    except:
        print('Error in combining audio and video')
        errorlog.append('Error in combining audio and video for file: ' + file)
        continue

# Now the same but with 16kHz audio
wavfolder16 = os.path.join(curfolder, 'data', 'Data_processed', 'Data_trials', 'Audio')
wavfiles16 = glob.glob(os.path.join(wavfolder16, '*.wav'))

for file in wavfiles16:

    print('Now processing file '+file)

    filename = file.split(os.sep)[-1].split('.')[0]
    dyadIndex = filename.split('_')[0]   # this is dyad number
    partIndex = filename.split('_')[1]

    if 'rein' in filename:
        sessionIndex = dyadIndex + '_' + partIndex + '_rein'
    else:
        sessionIndex = dyadIndex + '_' + partIndex

    if sessionIndex not in with16:
        continue

    if '_rein_' in filename: 
        trialIndex = filename.split('_')[4] # this is trial number
        participant = filename.split('_')[8]
        word = filename.split('_')[9]
        modality = filename.split('_')[10]
        
    else:
        trialIndex = filename.split('_')[3] # this is trial number
        participant = filename.split('_')[7] # this is participant 0/1
        word = filename.split('_')[8] # this is the word
        modality = filename.split('_')[9] # this is the modality

    # Handle correction
    if 'c0' in file:
        correction = '_c0'
    elif 'c1' in file:
        correction = '_c1'
    elif 'c2' in file:
        correction = '_c2'
    else:
        correction = ''

    # Trial type
    if '_pr_' in file:
        trialtype = 'pr'
    else:
        trialtype = 'trial'

    # This is the location where the video will be saved
    output_path = os.path.abspath(os.path.join(audiovideo, f"{sessionIndex}_{trialtype}_{trialIndex}_{participant}_{word}_{modality}{correction}_final.avi"))
  
    if os.path.exists(output_path):
        print('Video already exists, skipping...')
        continue

    # Load in the audio
    print('Loading the audio')
    audio_path = file
    if not os.path.exists(audio_path):
        print(f"Audio file not found: {audio_path}")

    # Input the video with ffmpg
    input_audio = ffmpeg.input(audio_path)

    # Load in the video with matchich trialIndex and SessionIndex
    print('Loading the video')
    video_path = os.path.join(videofolder, f"{sessionIndex}_{trialtype}_{trialIndex}_{participant}_{word}_{modality}{correction}{cut}_video_raw.avi")

    if not os.path.exists(video_path):
        print(f"Video file not found: {video_path}")
        errorlog.append('Video file not found: ' + video_path)

    input_video = ffmpeg.input(video_path)
    
    # Combine the audio and video
    print('Combining audio and video')

    try:
        ffmpeg.concat(input_video, input_audio, v=1, a=1).output(
            output_path,
            vcodec='libx264',
            acodec='aac',
            video_bitrate='2M',         
            
            ).run(overwrite_output=True)
    except:
        print('Error in combining audio and video')
        errorlog.append('Error in combining audio and video for file: ' + file)
        continue
      
# Write the error log
errorlog = pd.DataFrame(errorlog, columns=['File'])
currendate = pd.to_datetime("today").strftime("%Y%m%d")
outfile = os.path.join(audiovideo, f'errorlog_combining_audio_video_{currendate}.csv')

errorlog.to_csv(outfile, index=False, header=False)
print('All done!')

Now we are ready to visually inspect all files