import numpy as np import matplotlib import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation import matplotlib.patches as patches import librosa import random import pandas as pd import threading import time from dataclasses import dataclass from scipy.signal import stft, butter, sosfiltfilt from scipy.fft import rfft, rfftfreq matplotlib.use('TkAgg') # Константы SOUND_SPEED = 343.2 # скорость звука (м/с) MIC_DISTANCE = 0.06 # расстояние между микрофонами (м) ROOM_WIDTH = 3.0 # ширина комнаты (м) ROOM_HEIGHT = 2.0 # высота комнаты (м) SAMPLE_RATE = 48000 # частота дискретизации (Гц) CHUNK = 16384 # размер буфера RMS_THRESHOLD = 0.05 # порог RMS для определения звука FREQUENCY_RANGE = (500, 2000) # диапазон частот для поиска доминирующей частоты (Гц) ANGLE_RESOLUTION = 1 # шаг угла для поиска (градусы) SILENCE_TIMEOUT = 2.0 # время в секундах для сохранения последнего угла PEAK_THRESHOLD = 0.5 # порог для проверки достоверности пика в MUSIC-спектре MOVE_INTERVAL = 0.5 # интервал перемещения источника звука (с) @dataclass class Microphone: """Класс для хранения информации о микрофоне""" x: float y: float @dataclass class SoundSource: """Класс для хранения информации об источнике звука""" x: float y: float class MusicDirectionFinder: def __init__(self, mic_distance: float, audio_file: str): """Инициализация определителя направления с использованием MUSIC""" self.mic_distance = mic_distance self.mic1 = Microphone(x=-mic_distance / 2, y=0.0) self.mic2 = Microphone(x=mic_distance / 2, y=0.0) self.sound_source = self._generate_random_sound_source() self.running = True self.current_angle = 0.0 self.sound_detected = False self.last_sound_time = 0 self.last_detected_angle = None self.show_arrow = False self.rms_left = 0.0 self.rms_right = 0.0 self.audio_data, self.sample_rate = self.load_audio(audio_file) self.audio_index = 0 self.noise_level = 0.001 # шум self.results = [] self.source_positions = [(self.sound_source.x, self.sound_source.y, 0.0)] self.last_move_time = time.time() print("MusicDirectionFinder инициализирован") def _generate_random_sound_source(self) -> SoundSource: """Генерация положения источника звука с последовательным проходом углов от -90 до 90 градусов с шагом 20 градусов""" if not hasattr(self, 'angle_ranges'): self.angle_ranges = list(range(-90, 91, 20)) self.current_range_idx = 0 self.angle_count = 0 start_angle = self.angle_ranges[self.current_range_idx] end_angle = start_angle + 20 if self.current_range_idx < len(self.angle_ranges) - 1 else 90 angle_deg = random.uniform(start_angle, end_angle) angle_rad = np.radians(angle_deg) distance = 1.5 # расстояние source_x = distance * np.sin(angle_rad) source_y = distance * np.cos(angle_rad) if abs(source_x) > ROOM_WIDTH / 2: scale = (ROOM_WIDTH / 2) / abs(source_x) source_x *= scale source_y *= scale if source_y > ROOM_HEIGHT: scale = ROOM_HEIGHT / source_y source_x *= scale source_y *= scale if source_y < 0: source_y = 0.0 source_x = 0.0 print( f"Генерация новой позиции: угол={angle_deg:.2f}°, x={source_x:.2f}, y={source_y:.2f}, расстояние={distance:.2f}, повторение {self.angle_count + 1}/10") self.angle_count += 1 if self.angle_count >= 10: self.angle_count = 0 self.current_range_idx = (self.current_range_idx + 1) % len(self.angle_ranges) return SoundSource(x=source_x, y=source_y) def load_audio(self, filename: str) -> tuple: """Загрузка аудиофайла""" try: audio_data, sample_rate = librosa.load(filename, sr=SAMPLE_RATE, mono=True) print(f"Аудиофайл загружен: {filename}, длина: {len(audio_data)}") return audio_data, sample_rate except Exception as e: raise ValueError(f"Ошибка загрузки аудиофайла: {e}") def calculate_distances(self) -> tuple: """Расчет расстояний от источника звука до микрофонов""" l1 = np.sqrt((self.sound_source.x - self.mic1.x) ** 2 + (self.sound_source.y - self.mic1.y) ** 2) l2 = np.sqrt((self.sound_source.x - self.mic2.x) ** 2 + (self.sound_source.y - self.mic2.y) ** 2) return l1, l2 def process_signals_with_delay(self, signal: np.ndarray) -> tuple: """Обработка сигналов с учетом временного сдвига и шума""" l1, l2 = self.calculate_distances() t1 = l1 / SOUND_SPEED t2 = l2 / SOUND_SPEED log = int(self.sample_rate * (t2 - t1)) if len(signal) < CHUNK: signal = np.pad(signal, (0, CHUNK - len(signal)), mode='constant') elif len(signal) > CHUNK: signal = signal[:CHUNK] if log >= 0: S1 = signal[abs(log):] if log != 0 else signal.copy() S2 = signal[:-abs(log)] if log != 0 else signal.copy() else: S1 = signal[:-abs(log)] if log != 0 else signal.copy() S2 = signal[abs(log):] min_length = min(len(S1), len(S2)) S1 = S1[:min_length] S2 = S2[:min_length] if len(S1) < CHUNK: S1 = np.pad(S1, (0, CHUNK - len(S1)), mode='constant') S2 = np.pad(S2, (0, CHUNK - len(S2)), mode='constant') noise1 = np.random.normal(0, self.noise_level, size=S1.shape) noise2 = np.random.normal(0, self.noise_level, size=S2.shape) S1 = S1 + noise1 S2 = S2 + noise2 return S1, S2, t1, t2 def capture_audio(self): """Эмуляция захвата аудиоданных из файла""" if self.audio_index + CHUNK >= len(self.audio_data): self.audio_index = 0 chunk = self.audio_data[self.audio_index:self.audio_index + CHUNK] self.audio_index += CHUNK signal1, signal2, t1, t2 = self.process_signals_with_delay(chunk) print(f"capture_audio: длина signal1={len(signal1)}, signal2={len(signal2)}") return signal1, signal2, t1, t2 def calculate_rms(self, signal: np.ndarray) -> float: """Вычисление RMS сигнала""" return np.sqrt(np.mean(signal ** 2)) def find_dominant_frequency(self, signal: np.ndarray) -> float: """Определение доминирующей частоты сигнала через FFT""" if len(signal) != CHUNK: if len(signal) < CHUNK: signal = np.pad(signal, (0, CHUNK - len(signal)), mode='constant') else: signal = signal[:CHUNK] fft_signal = np.abs(rfft(signal)) freqs = rfftfreq(CHUNK, 1 / SAMPLE_RATE) if len(fft_signal) != len(freqs): raise ValueError(f"Несоответствие размеров fft_signal ({len(fft_signal)}) и freqs ({len(freqs)})") mask = (freqs >= FREQUENCY_RANGE[0]) & (freqs <= FREQUENCY_RANGE[1]) fft_signal = fft_signal[mask] freqs = freqs[mask] if len(fft_signal) == 0: print("find_dominant_frequency: нет частот в диапазоне, возвращаем FREQUENCY_RANGE[0]") return FREQUENCY_RANGE[0] dominant_freq = freqs[np.argmax(fft_signal)] print(f"find_dominant_frequency: доминирующая частота = {dominant_freq:.1f} Hz") return dominant_freq def bandpass_filter(self, signal: np.ndarray, freq: float) -> np.ndarray: """Применение полосового фильтра к сигналу""" bandwidth = 100 lowcut = max(FREQUENCY_RANGE[0], freq - bandwidth / 2) highcut = min(FREQUENCY_RANGE[1], freq + bandwidth / 2) sos = butter(4, [lowcut, highcut], btype='band', fs=SAMPLE_RATE, output='sos') filtered_signal = sosfiltfilt(sos, signal) return filtered_signal def steering_vector(self, angle: float, freq: float) -> np.ndarray: """Вычисление steering vector для заданного угла и частоты""" tau = (self.mic_distance * np.sin(np.radians(angle))) / SOUND_SPEED phase_shift = 2 * np.pi * freq * tau return np.array([1, np.exp(1j * phase_shift)]) # Исправлено: -1j на 1j def music_spectrum(self, signal1: np.ndarray, signal2: np.ndarray, freq: float) -> tuple: """Вычисление спектра MUSIC для оценки DoA""" f, t, Zxx1 = stft(signal1, fs=SAMPLE_RATE, nperseg=CHUNK, noverlap=int(CHUNK * 3 / 4)) _, _, Zxx2 = stft(signal2, fs=SAMPLE_RATE, nperseg=CHUNK, noverlap=int(CHUNK * 3 / 4)) freq_idx = np.argmin(np.abs(f - freq)) X = np.vstack((Zxx1[freq_idx, :], Zxx2[freq_idx, :])) R = np.dot(X, X.conj().T) / X.shape[1] R += np.eye(R.shape[0]) * 1e-6 eigenvalues, eigenvectors = np.linalg.eigh(R) idx = np.argsort(eigenvalues)[::-1] eigenvectors = eigenvectors[:, idx] En = eigenvectors[:, -1:] angles = np.arange(-90, 91, ANGLE_RESOLUTION) music_spectrum = np.zeros(len(angles)) for i, angle in enumerate(angles): a = self.steering_vector(angle, freq) music_spectrum[i] = 1 / np.abs(np.dot(a.conj().T, np.dot(En @ En.conj().T, a))) music_spectrum = music_spectrum / np.max(music_spectrum) max_idx = np.argmax(music_spectrum) print(f"music_spectrum: максимум спектра при угле {angles[max_idx]:.1f}°, значение {music_spectrum[max_idx]:.4f}") if music_spectrum[max_idx] < PEAK_THRESHOLD: print("music_spectrum: пик ниже порога, возвращаем 0.0") return 0.0, angles, music_spectrum estimated_angle = angles[max_idx] print(f"music_spectrum: расчётный угол = {estimated_angle:.1f}°") return estimated_angle, angles, music_spectrum def run(self): """Обработка аудио в реальном времени""" print("Запуск потока run") while self.running: try: current_time = time.time() if current_time - self.last_move_time >= MOVE_INTERVAL: self.sound_source = self._generate_random_sound_source() self.source_positions.append((self.sound_source.x, self.sound_source.y, current_time)) print(f"Источник звука перемещен в: x={self.sound_source.x:.2f}, y={self.sound_source.y:.2f}") self.last_move_time = current_time left, right, t1, t2 = self.capture_audio() self.rms_left = self.calculate_rms(left) self.rms_right = self.calculate_rms(right) signal_diff = np.mean(np.abs(left - right)) new_sound_detected = (self.rms_left > RMS_THRESHOLD) or (self.rms_right > RMS_THRESHOLD) if new_sound_detected: print(f"RMS left: {self.rms_left:.4f}, RMS right: {self.rms_right:.4f}") print(f"Signal difference: {signal_diff:.4f}") freq = self.find_dominant_frequency(left) left_filtered = self.bandpass_filter(left, freq) right_filtered = self.bandpass_filter(right, freq) angle, _, _ = self.music_spectrum(left_filtered, right_filtered, freq) self.current_angle = angle self.last_detected_angle = angle self.last_sound_time = current_time self.sound_detected = True self.show_arrow = True true_dx = self.sound_source.x true_dy = self.sound_source.y true_angle = np.arctan2(true_dx, true_dy) * 180 / np.pi self.results.append({ 't1 (ms)': t1 * 1000, 't2 (ms)': t2 * 1000, 'Time Delay (ms)': (t2 - t1) * 1000, 'Detected Angle (°)': angle, 'True Angle (°)': true_angle, 'Source X': self.sound_source.x, 'Source Y': self.sound_source.y }) else: if self.last_detected_angle is not None and current_time - self.last_sound_time < SILENCE_TIMEOUT: self.current_angle = self.last_detected_angle self.sound_detected = False self.show_arrow = True else: self.sound_detected = False self.show_arrow = False time.sleep(CHUNK / self.sample_rate) except Exception as e: print(f"Ошибка в run: {e}") continue print("Поток run завершён") def get_coordinates_dataframe(self): """Создание датафрейма с координатами микрофонов и всех позиций источника""" data = { 'Object': ['Mic1', 'Mic2'] + [f'SoundSource_{i}' for i in range(len(self.source_positions))], 'X': [self.mic1.x, self.mic2.x] + [pos[0] for pos in self.source_positions], 'Y': [self.mic1.y, self.mic2.y] + [pos[1] for pos in self.source_positions], 'Time': [0.0, 0.0] + [pos[2] for pos in self.source_positions] } return pd.DataFrame(data) def get_results_dataframe(self): """Создание датафрейма с результатами""" return pd.DataFrame(self.results) def main(): try: finder = MusicDirectionFinder(MIC_DISTANCE, "my_recording1.wav") print("Аудиофайл загружен, частота дискретизации:", finder.sample_rate) thread = threading.Thread(target=finder.run) thread.start() print("Поток run запущен") fig, ax = plt.subplots(figsize=(10, 8)) ax.set_xlim(-ROOM_WIDTH / 2, ROOM_WIDTH / 2) ax.set_ylim(-0.5, ROOM_HEIGHT) ax.set_title("Определение направления на источник звука (MUSIC)") ax.set_xlabel("X (м)") ax.set_ylabel("Y (м)") ax.grid(True) print("График инициализирован") ax.plot(finder.mic1.x, finder.mic1.y, 'bs', markersize=10, label='Микрофон 1') ax.plot(finder.mic2.x, finder.mic2.y, 'bs', markersize=10, label='Микрофон 2') source_plot, = ax.plot(finder.sound_source.x, finder.sound_source.y, 'ro', markersize=10, label='Источник звука') arrow_length = min(ROOM_WIDTH, ROOM_HEIGHT) / 3 arrow = ax.arrow(0, 0, 0, 0, head_width=0.05, head_length=0.1, fc='r', ec='r', label='Расчетное направление') arrow_true = ax.arrow(0, 0, 0, 0, head_width=0.05, head_length=0.1, fc='g', ec='g', linestyle=':', label='Истинное направление') sound_bar = ax.axhline(y=ROOM_HEIGHT - 0.2, color='green', linewidth=20, visible=False) angle_text = ax.text(0, ROOM_HEIGHT - 0.3, "", ha='center', va='center', fontsize=12) left_indicator = patches.Rectangle((finder.mic1.x - 0.02, -0.1), 0.04, 0.05, facecolor='gray') right_indicator = patches.Rectangle((finder.mic2.x - 0.02, -0.1), 0.04, 0.05, facecolor='gray') ax.add_patch(left_indicator) ax.add_patch(right_indicator) def update(frame): try: source_plot.set_data([finder.sound_source.x], [finder.sound_source.y]) dx = finder.sound_source.x dy = finder.sound_source.y true_angle_rad = np.arctan2(dx, dy) true_end_x = arrow_length * np.sin(true_angle_rad) true_end_y = arrow_length * np.cos(true_angle_rad) arrow_true.set_data(x=0, y=0, dx=true_end_x, dy=true_end_y) if finder.show_arrow: angle = finder.current_angle calc_angle_rad = np.radians(angle) calc_end_x = arrow_length * np.sin(calc_angle_rad) calc_end_y = arrow_length * np.cos(calc_angle_rad) arrow.set_data(x=0, y=0, dx=calc_end_x, dy=calc_end_y) angle_text.set_text(f"Расчетный угол: {angle:.1f}°\nИстинный угол: {np.degrees(true_angle_rad):.1f}°") if finder.sound_detected: sound_bar.set_visible(True) ax.set_title("Активное обнаружение звука (MUSIC)") else: sound_bar.set_visible(False) ax.set_title("Последнее зафиксированное направление (MUSIC)") else: arrow.set_data(x=0, y=0, dx=0, dy=0) angle_text.set_text("") sound_bar.set_visible(False) ax.set_title("Звук не обнаружен") left_indicator.set_facecolor('green' if finder.rms_left > RMS_THRESHOLD else 'gray') right_indicator.set_facecolor('green' if finder.rms_right > RMS_THRESHOLD else 'gray') return [source_plot, arrow, arrow_true, sound_bar, left_indicator, right_indicator, angle_text] except Exception as e: print(f"Ошибка в update: {e}") return [] ani = FuncAnimation(fig, update, frames=None, interval=100, blit=True, cache_frame_data=False) print("FuncAnimation инициализирован") plt.legend(loc='upper right') print("Вызов plt.show()") plt.show() print("plt.show() выполнен") finder.running = False thread.join() print("Поток run остановлен") coords_df = finder.get_coordinates_dataframe() results_df = finder.get_results_dataframe() print("\nКоординаты микрофонов и всех позиций источника звука:") print(coords_df.to_string(index=False)) print("\nРезультаты вычислений:") print(results_df.to_string(index=False)) coords_df.to_csv('coordinates_music.csv', index=False) results_df.to_csv('results_music.csv', index=False) print("\nДанные сохранены в 'coordinates_music.csv' и 'results_music.csv'") except Exception as e: print(f"Ошибка в main: {e}") finder.running = False thread.join() if __name__ == "__main__": main()