import numpy as np import librosa from scipy.signal import butter, sawtooth, sosfilt import sys import math import threading import multiprocessing as mp from queue import Queue def bandpass_filter(data, lowcut, highcut, fsinput, order=4): sos = butter(order, [lowcut, highcut], btype='bandpass', output='sos', fs=fsinput) return sosfilt(sos, data) def lowpass_filter(data, cutoff, fsinput, order=4): sos = butter(order, cutoff, btype='lowpass', output='sos', fs=fsinput) return sosfilt(sos, data) def generate_optical_track(input_wav, output_bin, height, width, fps, num_threads=None,slewrate=20): # 自动检测CPU核心数 if num_threads is None: num_threads = mp.cpu_count() # 加载并预处理音频 data, fs = librosa.load(input_wav, sr=None, mono=False) # 多声道判断 ismultichan = False channel = 1 if len(data.shape) > 1: ismultichan = True channel = data.shape[0] print(f"光刃曝光图案生成部分将使用{num_threads}个线程") print(f"{channel}个声道，输出视频：高度{height}、宽度{width*channel}、帧率{fps}") # 需要对所有样本取最大值，对每个通道除以同一个系数，以保持立体声声像 # 数据归一化 data /= np.max(np.abs(data)) print("11") # 滤波处理 highcut = (height * fps) / 2 if highcut >= fs / 2: highcut = (fs / 2) - 1 filtered = bandpass_filter(data, 40, highcut, fs) filtered /= np.max(np.abs(filtered)) # 波形叠加自身包络，将交变波形转为正极性脉冲，并使光刃宽度最小化 envelope = lowpass_filter(np.abs(data), 80, fs) envelope /= np.max(np.abs(envelope)) filtered = filtered + envelope filtered -= np.min(filtered) / 3 filtered /= np.max(filtered) print("22") # 预计算三角波模板（每行） t_line = np.linspace(0, 1, width, endpoint=False) triangle_line = (sawtooth(2 * np.pi * t_line + np.pi, 0.5) + 1) / 2 # 创建线程安全队列 frame_queue = Queue() result_queue = Queue() # 单声道处理方法 if not ismultichan: print("1chan") # 生成三角波模板（每帧） triangle_frame_template = np.tile(triangle_line, height) # 计算总帧数和每帧采样数 samples_per_frame = height * width total_time = len(filtered) / fs n_frames = math.ceil(total_time * fps) # 创建并启动工作线程 threads = [] for _ in range(num_threads): t = threading.Thread( target=frame_processing_worker, args=(filtered, fs, fps, samples_per_frame, triangle_frame_template, frame_queue, result_queue,slewrate) ) t.daemon = True t.start() threads.append(t) # 将帧索引放入队列供工作线程处理 for frame_idx in range(n_frames): frame_queue.put(frame_idx) # 添加停止信号 for _ in range(num_threads): frame_queue.put(None) # 等待所有工作线程完成 frame_queue.join() # 收集处理结果并按帧顺序排序 frame_results = [None] * n_frames for _ in range(n_frames): frame_idx, binary_data = result_queue.get() frame_results[frame_idx] = binary_data # 按顺序写入输出文件 with open(output_bin, 'wb') as f_out: for binary_data in frame_results: if binary_data is not None: f_out.write(binary_data) # 多声道处理方法 else: print("mulchan") # 将输出图像视为每帧宽度width、高度height*channel，每行时分复用的信号 # 生成三角波模板（每帧） triangle_frame_template = np.tile(triangle_line, height * channel) # 计算总帧数和每帧采样数 samples_per_frame = height * width total_time = len(filtered[0]) / fs n_frames = math.ceil(total_time * fps) # 创建并启动工作线程 threads = [] for _ in range(num_threads): t = threading.Thread( target=multichannel_frame_processing_worker, args=(filtered, fs, fps, samples_per_frame, channel, triangle_frame_template, frame_queue, result_queue, height, width,slewrate) ) t.daemon = True t.start() threads.append(t) # 将帧索引放入队列供工作线程处理 for frame_idx in range(n_frames): print(frame_idx,"/",n_frames) frame_queue.put(frame_idx) # 添加停止信号 for _ in range(num_threads): frame_queue.put(None) # 等待所有工作线程完成 frame_queue.join() # 收集处理结果并按帧顺序排序 frame_results = [None] * n_frames for _ in range(n_frames): frame_idx, binary_data = result_queue.get() frame_results[frame_idx] = binary_data # 按顺序写入输出文件 with open(output_bin, 'wb') as f_out: for binary_data in frame_results: if binary_data is not None: f_out.write(binary_data) def frame_processing_worker(filtered, fs, fps, samples_per_frame, triangle_frame_template, frame_queue, result_queue,slewrate): print(result_queue) fs_fps_ratio = fs / fps while True: frame_idx = frame_queue.get() if frame_idx is None: frame_queue.task_done() break # 计算当前帧所有采样点在原信号中的位置 j = np.arange(samples_per_frame) pos = (frame_idx + j / samples_per_frame) * fs_fps_ratio # 初始化当前帧的重采样信号 frame_resampled = np.zeros(samples_per_frame) # 有效位置处理 (0 <= pos < len(filtered) - 1) mask_valid = (pos >= 0) & (pos < len(filtered) - 1) valid_pos = pos[mask_valid] if len(valid_pos) > 0: # 线性插值 idx0 = valid_pos.astype(int) idx1 = idx0 + 1 frac = valid_pos - idx0 frame_resampled[mask_valid] = (1 - frac) * filtered[idx0] + frac * filtered[idx1] # 使用预计算的三角波模板计算输出 # 要是去掉了三角波比较，直接输出frame_resampled，就成了变密式音轨 output_frame = slewrate * (frame_resampled - triangle_frame_template) # output_frame = frame_resampled output_frame = np.clip(output_frame, 0, 1) binary_frame = (output_frame * 255).astype(np.uint8) # 将结果放入结果队列 result_queue.put((frame_idx, binary_frame.tobytes())) frame_queue.task_done() def multichannel_frame_processing_worker(filtered, fs, fps, samples_per_frame, channel, triangle_frame_template, frame_queue, result_queue, height, width,slewrate): print(result_queue) fs_fps_ratio = fs / fps total_samples = samples_per_frame * channel while True: frame_idx = frame_queue.get() if frame_idx is None: frame_queue.task_done() break # 计算当前帧所有采样点在原信号中的位置（向量化） j = np.arange(samples_per_frame) pos = (frame_idx + j / samples_per_frame) * fs_fps_ratio # 初始化当前帧的重采样信号 frame_resampled = np.zeros((channel, samples_per_frame)) # 有效位置处理 (0 <= pos < len(filtered[0]) - 1) mask_valid = (pos >= 0) & (pos < len(filtered[0]) - 1) valid_pos = pos[mask_valid] if len(valid_pos) > 0: # 向量化线性插值 idx0 = valid_pos.astype(int) idx1 = idx0 + 1 frac = valid_pos - idx0 # 对所有声道进行向量化插值 for ch in range(channel): # 使用NumPy高级索引进行向量化操作 chan_data = filtered[ch] frame_resampled[ch, mask_valid] = ( (1 - frac) * chan_data[idx0] + frac * chan_data[idx1] ) # 处理无效位置 frame_resampled[:, ~mask_valid] = 0 # 高效交织：使用reshape和transpose代替手动索引计算 # 1. 将数据重塑为 (channel, height, width) reshaped = frame_resampled.reshape((channel, height, width)) # 2. 转置为 (height, channel, width) transposed = reshaped.transpose(1, 0, 2) # 3. 重塑为一维数组 (height * channel * width) frame_resampled_interleaved = transposed.reshape(-1) # 使用预计算的三角波模板计算输出 # 要是去掉了三角波比较，直接输出frame_resampled，就成了变密式音轨 output_frame = slewrate * (frame_resampled_interleaved - triangle_frame_template) # output_frame = frame_resampled_interleaved output_frame = np.clip(output_frame, 0, 1) binary_frame = (output_frame * 255).astype(np.uint8) # 将结果放入结果队列 result_queue.put((frame_idx, binary_frame.tobytes())) frame_queue.task_done() if __name__ == '__main__': if len(sys.argv) < 6: print("双边变积式光学音轨仿真程序") print("使用方法：python opttkgen.py 输入文件名（使用librosa库，不局限于wav格式） 输出文件名（RAWVIDEO，Gray（uint8）） 每帧图像高度 每声道宽度 视频帧率 [线程数]") print("线程数不输入时默认为CPU线程数") sys.exit(1) input_wav = sys.argv[1] output_bin = sys.argv[2] height = int(sys.argv[3]) width = int(sys.argv[4]) fps = int(sys.argv[5]) # 默认使用使用CPU核心数，如果提供了第6个参数则使用它 num_threads = None if len(sys.argv) > 6: num_threads = int(sys.argv[6]) # 模仿真实比较器的摆率，数值越小，音轨图案的边缘越平滑、模糊。建议每声道宽度为64时设置摆率为20，音轨越宽，设置的值应越大。 slewrate=10 print("双边变积式光学音轨仿真程序") print("本程序由TCJY（核心算法、程序排错）、DeepSeek（内存占用优化、多声道、多线程部分）、ChatGPT（程序框架、单声道部分）合作编写") generate_optical_track(input_wav, output_bin, height, width, fps, num_threads,slewrate) print(f"输出已写入 {output_bin}。")