Полное руководство по устранению проблем синхронизации измерений ESP32 INMP441

import numpy as np
import serial
import time
import sounddevice as sd
from scipy.signal import chirp, correlate, windows

# === Улучшенная конфигурация ===
FS = 48000
DURATION = 0.0625
N_SAMPLES = int(FS * DURATION)
BUFFER_EXTRA = 2000  # Дополнительные выборки для гарантии захвата
SERIAL_PORT = 'COM8'
BAUD_RATE = 115200  # Стандартная скорость передачи
SYNC_DELAY = 0.1  # Задержка между готовностью ESP32 и разверткой

def generate_sweep():
    t = np.linspace(0, DURATION, N_SAMPLES, endpoint=False)
    f0 = 300
    f1 = 15000
    beta = np.log(f1 / f0) / DURATION
    phase = 2 * np.pi * f0 * (np.exp(beta * t) - 1) / beta
    sweep = np.sin(phase) * windows.hann(N_SAMPLES)
    return sweep

# === Улучшенная синхронизация ===
ser = serial.Serial(SERIAL_PORT, BAUD_RATE, timeout=5)
print("Ожидание сигнала READY от ESP32...")
while True:
    line = ser.readline().decode('utf-8').strip()
    if line == "READY":
        print("ESP32 готов. Воспроизведение развертки через 100мс...")
        time.sleep(SYNC_DELAY)  # Точная временная задержка
        break

# === Воспроизведение развертки с мониторингом ===
sweep = generate_sweep()
print(f"Воспроизведение развертки от {300}Гц до {15000}Гц...")
sd.play(sweep, samplerate=FS)
sd.wait()

# === Многоканальная запись ===
def record_with_loopback():
    print("Запись аудио с петлевой обратной связью...")
    recording = sd.rec(int((DURATION + 0.1) * FS), samplerate=FS, channels=2, blocking=True)
    reference_signal = recording[:, 0]  # Канал воспроизведения
    mic_signal = recording[:, 1]       # Входной канал ESP32
    return reference_signal, mic_signal

# Используйте это вместо последовательной связи для синхронизации

// Код Arduino для ESP32 - внутренняя генерация развертки
#include <driver/i2s.h>

#define I2S_WS 22
#define I2S_SD 12
#define I2S_SCK 21
#define I2S_PORT I2S_NUM_0

void setup() {
    Serial.begin(115200);
    i2s_config_t i2s_config = {
        .mode = (i2s_mode_t)(I2S_MODE_MASTER | I2S_MODE_RX),
        .sample_rate = 48000,
        .bits_per_sample = I2S_BITS_PER_SAMPLE_32BIT,
        .channel_format = I2S_CHANNEL_FMT_ONLY_LEFT,
        .communication_format = I2S_COMM_FORMAT_STAND_I2S,
        .intr_alloc_flags = ESP_INTR_FLAG_LEVEL1,
        .dma_buf_count = 8,
        .dma_buf_len = 64,
        .use_apll = true,  // Использовать APLL для лучшей точности тактовой частоты
        .tx_desc_auto_clear = false,
        .fixed_mclk = 0
    };
    
    i2s_pin_config_t pin_config = {
        .bck_io_num = I2S_SCK,
        .ws_io_num = I2S_WS,
        .data_out_num = I2S_PIN_NO_CHANGE,
        .data_in_num = I2S_SD
    };
    
    i2s_driver_install(I2S_PORT, &i2s_config, 0, NULL);
    i2s_set_pin(I2S_PORT, &pin_config);
}

void loop() {
    // Внутренняя генерация и запись развертки
    generate_and_measure_sweep();
}

# === Улучшенная взаимная корреляция ===
def align_signals(reference, recorded, max_lag_samples=1000):
    # Нормализация сигналов для лучшей корреляции
    ref_norm = (reference - np.mean(reference)) / np.std(reference)
    rec_norm = (recorded - np.mean(recorded)) / np.std(recorded)
    
    # Вычисление корреляции с несколькими методами
    corr = correlate(rec_norm, ref_norm, mode='full', method='auto')
    
    # Поиск пика с точностью до подвыборки с помощью параболической интерполяции
    peak_idx = np.argmax(corr)
    if peak_idx > 0 and peak_idx < len(corr) - 1:
        y1, y2, y3 = corr[peak_idx-1], corr[peak_idx], corr[peak_idx+1]
        peak_idx += 0.5 * (y1 - y3) / (2 * y2 - y1 - y3)
    
    lag_samples = int(peak_idx) - len(reference) + 1
    return lag_samples

def calculate_snr(signal, noise_samples=500):
    noise_region = signal[:noise_samples]
    signal_region = signal[noise_samples:2*noise_samples]
    noise_power = np.mean(noise_region ** 2)
    signal_power = np.mean(signal_region ** 2)
    snr_db = 10 * np.log10(signal_power / (noise_power + 1e-12))
    return snr_db

def detect_sweep_presence(signal, threshold=0.1):
    # Проверка, содержит ли сигнал ожидаемые характеристики развертки
    fft_signal = np.fft.rfft(signal)
    freqs = np.fft.rfftfreq(len(signal), 1/FS)
    
    # Поиск энергии в ожидаемом частотном диапазоне
    mask = (freqs >= 300) & (freqs <= 15000)
    energy_in_band = np.sum(np.abs(fft_signal[mask]) ** 2)
    total_energy = np.sum(np.abs(fft_signal) ** 2)
    
    return (energy_in_band / total_energy) > threshold

# Измерение фактической частоты дискретизации из записанных данных
def measure_actual_sample_rate(signal, expected_freq=1000):
    # Генерация тестового тона известной частоты
    test_tone = np.sin(2 * np.pi * expected_freq * np.arange(len(signal)) / FS)
    
    # Взаимная корреляция для определения фактической частоты
    corr = correlate(signal, test_tone, mode='same')
    peak_idx = np.argmax(corr)
    
    # Расчет фактической частоты дискретизации на основе положения пика
    actual_rate = FS * expected_freq / (expected_freq + (peak_idx - len(signal)//2) * expected_freq / len(signal))
    return actual_rate

def apply_windowed_fft(signal, window_type='hann', overlap=0.5):
    window = windows.get_window(window_type, len(signal))
    windowed_signal = signal * window
    
    # Применение коэффициента коррекции окна
    window_correction = np.sum(window ** 2) / len(window)
    corrected_signal = windowed_signal / window_correction
    
    return np.fft.rfft(corrected_signal)

def remove_dc_offset_enhanced(signal, method='highpass'):
    if method == 'highpass':
        # Простой фильтр высоких частот
        alpha = 0.99
        filtered = np.zeros_like(signal)
        filtered[0] = signal[0]
        for i in range(1, len(signal)):
            filtered[i] = alpha * filtered[i-1] + (1-alpha) * (signal[i] - signal[i-1])
        return filtered
    elif method == 'polynomial':
        # Подгонка и удаление полиномиальной базовой линии
        x = np.arange(len(signal))
        coeffs = np.polyfit(x, signal, deg=3)
        baseline = np.polyval(coeffs, x)
        return signal - baseline

// Улучшенная конфигурация I2S с APLL
i2s_config_t i2s_config = {
    .mode = (i2s_mode_t)(I2S_MODE_MASTER | I2S_MODE_RX),
    .sample_rate = 48000,
    .bits_per_sample = I2S_BITS_PER_SAMPLE_32BIT,
    .channel_format = I2S_CHANNEL_FMT_ONLY_LEFT,
    .communication_format = I2S_COMM_FORMAT_STAND_I2S,
    .intr_alloc_flags = ESP_INTR_FLAG_LEVEL1,
    .dma_buf_count = 8,
    .dma_buf_len = 64,
    .use_apll = true,  // Включить APLL для лучшей точности тактовой частоты
    .tx_desc_auto_clear = false,
    .fixed_mclk = 0
};

# Динамическое изменение размера буфера на основе длительности развертки
BUFFER_DURATION = 0.2  # 200мс для гарантии захвата
BUFFER_SAMPLES = int(FS * BUFFER_DURATION)

def process_inmp441_data(raw_data):
    # Преобразование 32-битного знакового в 24-битный эффективно
    # Сдвиг вправо на 8 бит для правильного выравнивания (может потребовать коррекции в зависимости от вашего конкретного устройства)
    processed = np.array(raw_data, dtype=np.int32) >> 8
    
    # Преобразование в float и нормализация
    processed = processed.astype(np.float32) / (2**23 - 1)
    
    # Применение коррекции частотной характеристики INMP441, если известно
    # (Это потребовало бы калибровочных данных)
    return processed

import numpy as np
import matplotlib.pyplot as plt
import serial
import time
import sounddevice as sd
from scipy.signal import chirp, correlate, windows
from scipy import signal as scipy_signal

class FrequencyResponseAnalyzer:
    def __init__(self, fs=48000, duration=0.0625, extra_buffer=0.1):
        self.fs = fs
        self.duration = duration
        self.extra_buffer = extra_buffer
        self.n_sweep_samples = int(fs * duration)
        self.n_buffer_samples = int(fs * (duration + extra_buffer))
        
    def generate_sweep(self, f0=300, f1=15000):
        """Генерация логарифмической частотной развертки с окном Ханна"""
        t = np.linspace(0, self.duration, self.n_sweep_samples, endpoint=False)
        beta = np.log(f1 / f0) / self.duration
        phase = 2 * np.pi * f0 * (np.exp(beta * t) - 1) / beta
        sweep = np.sin(phase) * windows.hann(self.n_sweep_samples)
        return sweep
    
    def record_with_esp32(self, serial_port='COM8', baud_rate=115200):
        """Улучшенная запись с ESP32 с лучшим таймингом"""
        ser = serial.Serial(serial_port, baud_rate, timeout=5)
        
        # Ожидание готовности
        print("Ожидание готовности ESP32...")
        while True:
            line = ser.readline().decode('utf-8').strip()
            if line == "READY":
                break
        
        # Воспроизведение развертки с точным таймингом
        sweep = self.generate_sweep()
        print("Воспроизведение развертки...")
        sd.play(sweep, samplerate=self.fs)
        sd.wait()
        
        # Запись данных
        print("Запись ответа...")
        raw_data = []
        while len(raw_data) < self.n_buffer_samples:
            if ser.in_waiting:
                line = ser.readline().decode('utf-8').strip()
                try:
                    raw_data.append(float(line))
                except ValueError:
                    continue
        
        ser.close()
        return np.array(raw_data)
    
    def align_signals(self, reference, recorded):
        """Улучшенное выравнивание сигналов с точностью до подвыборки"""
        # Нормализация сигналов
        ref_norm = (reference - np.mean(reference)) / (np.std(reference) + 1e-12)
        rec_norm = (recorded - np.mean(recorded)) / (np.std(recorded) + 1e-12)
        
        # Вычисление корреляции
        corr = correlate(rec_norm, ref_norm, mode='full', method='auto')
        
        # Поиск пика с точностью до подвыборки
        peak_idx = np.argmax(corr)
        if 0 < peak_idx < len(corr) - 1:
            y1, y2, y3 = corr[peak_idx-1], corr[peak_idx], corr[peak_idx+1]
            peak_idx += 0.5 * (y1 - y3) / (2 * y2 - y1 - y3)
        
        lag_samples = int(peak_idx) - len(reference) + 1
        return lag_samples
    
    def validate_measurement(self, signal, min_snr_db=20):
        """Проверка качества измерения"""
        # Расчет SNR
        noise_region = signal[:500]
        signal_region = signal[500:1000]
        noise_power = np.mean(noise_region ** 2)
        signal_power = np.mean(signal_region ** 2)
        snr_db = 10 * np.log10(signal_power / (noise_power + 1e-12))
        
        # Проверка наличия развертки
        fft_signal = np.fft.rfft(signal)
        freqs = np.fft.rfftfreq(len(signal), 1/self.fs)
        mask = (freqs >= 300) & (freqs <= 15000)
        energy_in_band = np.sum(np.abs(fft_signal[mask]) ** 2)
        total_energy = np.sum(np.abs(fft_signal) ** 2)
        sweep_present = (energy_in_band / total_energy) > 0.1
        
        return {
            'snr_db': snr_db,
            'sweep_present': sweep_present,
            'valid': snr_db >= min_snr_db and sweep_present
        }
    
    def calculate_frequency_response(self, reference, recorded, nfft=8192):
        """Расчет частотной характеристики с улучшенной обработкой"""
        # Удаление смещения постоянной составляющей
        reference = scipy_signal.detrend(reference)
        recorded = scipy_signal.detrend(recorded)
        
        # Применение оконного преобразования
        window = windows.hann(len(reference))
        reference_windowed = reference * window
        recorded_windowed = recorded[:len(reference)] * window
        
        # Вычисление FFT
        S = np.fft.rfft(reference_windowed, nfft)
        M = np.fft.rfft(recorded_windowed, nfft)
        
        # Расчет передаточной функции
        H = M / (S + 1e-12)
        
        # Преобразование в децибелы
        freqs = np.fft.rfftfreq(nfft, 1/self.fs)
        mag_db = 20 * np.log10(np.abs(H) + 1e-12)
        phase_deg = np.angle(H, deg=True)
        
        return freqs, mag_db, phase_deg
    
    def run_measurement(self):
        """Полный рабочий процесс измерения"""
        # Запись данных
        recorded_data = self.record_with_esp32()
        
        # Генерация эталонной развертки
        reference_sweep = self.generate_sweep()
        
        # Выравнивание сигналов
        lag = self.align_signals(reference_sweep, recorded_data)
        aligned_recorded = recorded_data[lag:lag + len(reference_sweep)]
        
        # Проверка измерения
        validation = self.validate_measurement(aligned_recorded)
        if not validation['valid']:
            print(f"⚠️ Проверка измерения не пройдена: SNR {validation['snr_db']:.1f}дБ, Развертка присутствует: {validation['sweep_present']}")
            return None
        
        # Расчет частотной характеристики
        freqs, mag_db, phase_deg = self.calculate_frequency_response(reference_sweep, aligned_recorded)
        
        # Построение графиков
        self.plot_results(freqs, mag_db, phase_deg)
        
        return freqs, mag_db, phase_deg
    
    def plot_results(self, freqs, mag_db, phase_deg):
        """Построение графика частотной характеристики с информацией о проверке"""
        plt.figure(figsize=(12, 8))
        
        # График амплитуды
        plt.subplot(2, 1, 1)
        plt.semilogx(freqs, mag_db, 'b-', linewidth=2)
        plt.axhline(84.0, color='gray', linestyle='--', alpha=0.7, label='Эталонный УЗД (84 дБ)')
        plt.axhline(0.0, color='red', linestyle='--', alpha=0.7, label='0 дБ эталон')
        plt.ylabel('Амплитуда [дБ УЗД]')
        plt.title('Измерение частотной характеристики INMP441')
        plt.grid(True, which='both', alpha=0.3)
        plt.legend()
        plt.xlim([20, 20000])
        
        # График фазы
        plt.subplot(2, 1, 2)
        plt.semilogx(freqs, phase_deg, 'g-', linewidth=2)
        plt.xlabel('Частота [Гц]')
        plt.ylabel('Фаза [градусы]')
        plt.grid(True, which='both', alpha=0.3)
        plt.xlim([20, 20000])
        
        plt.tight_layout()
        plt.show()

# Пример использования
if __name__ == "__main__":
    analyzer = FrequencyResponseAnalyzer()
    results = analyzer.run_measurement()
    if results:
        freqs, mag_db, phase_deg = results
        print("Измерение успешно завершено!")
        # Здесь можно выполнить дополнительный анализ