48.0KHz 的音频降采样为44.1KHz，明显不太容易，因为它并不能通过简单的抽样来实现。以标题为引，本文下述内容将全面的讲解重采样技术并提供实现代码。并且比MATLAB和Python/scipy 提供的resample函数效率更高，便于进行嵌入式开发。可以先收藏一下，以后没准用的到。
重采样有下述5种情况：
① 整倍数降采样；
② 整数倍升采样；
③ 分数比例重采样
④ 任意比例重采样
⑤ 角度域重采样
一、降采样：
在振动数据和音频数据降采样时，首先是不能直接抽取，因为会产生混叠现象。下图为320Hz 正弦数据，采样率1000Hz,直接进行5倍抽取至200Hz，从时域观察看出现了一个不应该存在的周期信号。看下面的频谱图，该频率为80Hz。降采样产生的无中生有的“混叠”信号，可能会让人产生误判。所以在抽取数据前先滤波。
根据奈奎斯特定律：200Hz采样率内不应该有大于100Hz的信号，因此降采样前采用低通为100Hz的低通滤波器对1000Hz数据进行滤波。
本文利用Sinc函数设计低通滤波器
num_taps = 2048 # sinc 波形长度，越长越好。if l <= num_taps: # 如果信号长度小于滤波器长度，则将num_taps设置为小于信号长度的最大的可以被2整除的数num_taps = l-1 if l % 2 == 1 else lt = np.arange(-num_taps // 2, num_taps // 2) / fs # fs 原始信号采样率，如上例中的1000Hzsinc_filter = np.sinc(2 * cutoff_freq * t) # cutoff_freq 截断频率B，如上例中的100Hz。因为归一化所以，没有乘幅值系数2Bsinc_filter *= get_window('hamming', num_taps) # 应用汉明窗截断sinc函数sinc_filter /= np.sum(sinc_filter) # 归一化result = np.convolve(a, sinc_filter, mode='same') #卷积
note: “这个滤波器的设计原理将另起一篇文章讲解”
公式参考："https://en.wikipedia.org/wiki/Sinc_filter”
二、升采样：
升采样是不是就是直接插值？答案是否定的。用差值代替重采样会引入噪音和误差：例：对长度0.02s,采样率为1000Hz，幅值为1，频率为320Hz的正弦信号进行重采样和插值到2000Hz.
note: “上述绘图没有加窗所以有一些频谱泄露，幅值小于1，频谱泄露不是本文重点”
插值后：从时域来看虽然波形与原始波形一致，但是在频谱上却使幅值降低明显，同时在高频处引入噪音。音频质量下降。
重采样：从时域来看，信号有一些因为滤波造成的延迟（当然这个延迟是容易弥补的），但很好的绘制了该正弦波形。在频率上没有引入噪音，音频质量提高。
重采样技术的实现与线性插值不同，是通过先插0后滤波来实现的。
低通滤波器截至带宽为：bw = 原始数据采样fs/2。
同样采用Sinc函数设计低通滤波器：
num_taps = 2048 # sinc 波形长度，越长越好。t = np.arange(-num_taps // 2, num_taps // 2) / fs_new # 升采样的目标采样率 ，上例2000hzsinc_filter = np.sinc(2 * cutoff_freq * t) # cutoff_freq = bw = 原始数据采样fs/2。sinc_filter *= get_window('hamming', num_taps) # 应用汉明窗sinc_filter /= np.sum(sinc_filter) # 归一化x_new = factor * np.convolve(x_new, sinc_filter, mode='same') # 卷积滤波
上述讲述降采样与升采样的基本过程，但是因为使用的插0与抽取数据，所以上述方法仅适用于整数倍的采样率处理。
如何和来完成标题中“48.0KHz 的音频重采样为44.1KHz”的任务呢？
方法一：先将48000Hz,升采样147倍至两个采样率的最小公倍数7056000，然后降采样160倍至44100Hz。这种方法是MATLAB采用的方法，即分数比例重采样，有时计算开销比较大。。。。。
方法二：傅里叶变换重采样，先将原始数据进行fft变换、在频域中进行数据处理然后进行傅里叶逆变换。Python/scipy 提供的函数resample即采用这种方法。稳定的计算开销较大。。。。。。
方法三：插值法
流程如下：
红圈为原始数据，采样率为48000Hz；首先进行15倍升采样，获取蓝色圆点数据，采样率48000Hz*15；然后进行线性插值至44100Hz采样点数据，如上图绿色数据点。完整代码再文章末尾。
factor = 15 # 15倍插0x_new = np.zeros(len(x) * factor)# 插零for i in range(len(x)):x_new[i * factor] = x[i]# 滤波cutoff_freq = fs_new * cutoff_freq_pre # 此处滤波器截止频率为目标采样率的一半 (Hz) cutoff_freq_pre默认滤波0.5t = np.arange(-num_taps // 2, num_taps // 2) / (fs * factor) # 目标升采样率(fs * upsample_factor)sinc_filter = np.sinc(2 * cutoff_freq * t)sinc_filter *= get_window('hamming', num_taps) # 应用汉明窗sinc_filter /= np.sum(sinc_filter)x_15 = factor * np.convolve(x_new, sinc_filter, mode='same') # 差入了factor个0， 需要补偿回来t_15 = np.arange(0, len(x_new), 1) / (fs * factor)t_new = np.arange(0, round(len(x_15)*fs_new/(factor*fs)), 1) / fs_newx_new = np.interp(t_new, t_15, x_15)
补充：上述过程在进行15倍升采样时滤波器截至频率为44100Hz/2;
线性插值时可先插值到采样率44100Hz的倍数，再进一步降采样；
误差分析：
SDR，全称是Signal to Distortion Ratio，中文可以翻译为信噪比。在信号处理中，信噪比是用来衡量信号质量的一个指标，它表示的是信号功率与噪声功率的比值。在SDR中，信号是指我们想要处理的原始数据，而噪声是指信号处理过程中引入的干扰和误差。信噪比越高，说明信号的质量越好，噪声的影响越小。SDR的单位是分贝（dB）
限于篇幅，本文对于角度域重采样先不涉及。
代码参考：下述代码支持任何频率的重采样。通时该功能已经集成到了本号（deepdaq）的开源振动噪音软件中：
import mathimport timeimport numpy as npfrom matplotlib import pyplot as pltfrom matplotlib.ticker import MultipleLocatorfrom scipy.signal import get_window# 升采样def _resample(x, fs, fs_new, cutoff_freq_pre, num_taps=2048):x_new = []t_new = []l = len(x) # 信号长度if l <= num_taps: # 如果信号长度小于滤波器长度，则将num_taps设置为小于信号长度的最大的可以被2整除的数num_taps = l - 1 if l % 2 == 1 else lif fs < fs_new: # 升采样if fs_new % fs == 0: # 整数倍升采样print(fs,fs_new,'ok')factor = fs_new // fsx_new = np.zeros(len(x) * factor)# 插零for i in range(len(x)):x_new[i*factor] = x[i]# 滤波cutoff_freq = fs * cutoff_freq_pre # 滤波器截止频率 (Hz) 默认滤波0.5# 对插0后的函数进行滤波，因为这个滤波器是先插值后滤波所以这个截至频率计算公式为 fs* cutoff_freq_pre/fs_new，所以t要除以fs_newt = np.arange(-num_taps // 2, num_taps // 2) / fs_newsinc_filter = np.sinc(2 * cutoff_freq * t)sinc_filter *= get_window('hamming', num_taps) # 应用汉明窗sinc_filter /= np.sum(sinc_filter)x_new = factor * np.convolve(x_new, sinc_filter, mode='same')t_new = np.arange(0, len(x_new), 1) / fs_newelse: # 非整数倍升采样# 分数倍升采样，计算最小公倍数，如果最小公倍数大于原始低样率的15倍，则为任意采样率比# 计算最大公约数gcd = math.gcd(fs, fs_new)# 计算最小公倍数lcm = (fs * fs_new) // gcdif lcm // fs <= 50: # 分数倍升采样，先将原始采样率升采样至最小公倍数，一般不应该超过50倍print((fs * fs_new)//math.gcd(fs, fs_new) // fs_new)# 先升采样至最大公倍数print(fs,lcm)x_new, t_new = _resample(x, fs, lcm, 0.5)# 在降采样x_new, t_new = _resample(x_new, lcm, fs_new, 0.5)else: # 如果重采样前后采样率接近,fs_new // fs == 1，则进行15倍升采样（提升精度）,插值成目标采样率。if fs_new // fs == 1:upsample_factor = 15x_15, t_15 = _resample(x, fs, upsample_factor*fs, 0.5) # 已经低通滤波了t_new = np.arange(0, round(len(x_15)*fs_new/(upsample_factor*fs)), 1) / fs_newx_new = np.interp(t_new, t_15, x_15)elif fs_new // fs > 1: #如果重采样前后采样率不接近,fs_new // fs > 1，先二倍升采样，然后进行15倍升采样（提升精度）,插值成目标采样率。print('run')x_2, t_2 = _resample(x, fs, 2 * fs, 0.5) # 此处进入递归x_new, t_new = _resample(x_2, 2 * fs, fs_new, 0.5)elif fs > fs_new: # 降采样if fs % fs_new == 0: # 整数倍降采样factor = fs // fs_new# 滤波cutoff_freq = fs_new * cutoff_freq_pre # 滤波器截止频率 (Hz) 默认滤波0.5# 对原始函数进行滤波，因为这个滤波器是先滤波后插值所以这个截至频率计算公式为 fs_new* cutoff_freq_pre/fs，所以t要除以fst = np.arange(-num_taps // 2, num_taps // 2) / fs # 对原始函数进行滤波sinc_filter = np.sinc(2 * cutoff_freq * t)sinc_filter *= get_window('hamming', num_taps) # 应用汉明窗sinc_filter /= np.sum(sinc_filter)x_new = np.convolve(x, sinc_filter, mode='same')# 抽样x_new = x_new[::factor]# x_new = x[::factor] # 直接抽取t_new = np.arange(0, len(x_new), 1) / fs_newelse: # 非整数倍降采样# 分数倍降采样，计算最小公倍数，如果最小公倍数大于原始低样率的15倍，则为任意采样率比# 计算最大公约数gcd = math.gcd(fs, fs_new)# 计算最小公倍数lcm = (fs * fs_new) // gcdif lcm // fs_new <= 50: # 分数倍升采样，先将原始采样率升采样至最小公倍数，一般不应该超过50倍print((fs * fs_new) // math.gcd(fs, fs_new) // fs_new)# 先升采样至最大公倍数print(fs, lcm)x_new, t_new = _resample(x, fs, lcm, 0.5)print(lcm, fs_new, 'run')# 在降采样x_new, t_new = _resample(x_new, lcm, fs_new, 0.5)else: #进行15倍升采样（提升精度）,插值成目标采样率。factor = 15 # 15倍插0x_new = np.zeros(len(x) * factor)# 插零for i in range(len(x)):x_new[i * factor] = x[i]# 滤波cutoff_freq = fs_new * cutoff_freq_pre # 此处滤波器截止频率为目标采样率的一半 (Hz) cutoff_freq_pre默认滤波0.5t = np.arange(-num_taps // 2, num_taps // 2) / (fs * factor) # 目标升采样率(fs * upsample_factor)sinc_filter = np.sinc(2 * cutoff_freq * t)sinc_filter *= get_window('hamming', num_taps) # 应用汉明窗sinc_filter /= np.sum(sinc_filter)x_15 = factor * np.convolve(x_new, sinc_filter, mode='same') # 差入了factor个0， 需要补偿回来t_15 = np.arange(0, len(x_new), 1) / (fs * factor)t_new = np.arange(0, round(len(x_15)*fs_new/(factor*fs)), 1) / fs_newx_new = np.interp(t_new, t_15, x_15)elif fs == fs_new: # 采样率不变x_new = xt_new = np.arange(0, len(x_new), 1) / fs_newreturn x_new, t_new# 绘制频率响应def plot_frequency_response(x, fs):# 计算频率响应n = len(x)frequencies = np.fft.rfftfreq(n, 1 / fs)fft_amplitude = 2*np.abs(np.fft.rfft(x))/nreturn frequencies, fft_amplitudeif __name__ == '__main__':fs = 1000 # 采样率duration = 0.02 # 信号持续时间，单位为秒n_samples = int(fs * duration)# 生成一个正弦信号t = np.linspace(0, duration, n_samples, endpoint=False)# 生成正弦信号x = 0*np.sin(2 * np.pi * 50 * t)+np.sin(2 * np.pi * 320 * t)# white_noise = np.random.normal(0, 1, n_samples)t_original = np.arange(0, n_samples, 1) / fs# 采样率加倍fs_new = 2000start_time = time.time()x_new, t_new = _resample(x, fs, fs_new, 0.5)end_time = time.time()print('resample time:', end_time - start_time)# 二倍插值t_interp = np.linspace(0, duration, n_samples * 2, endpoint=False)x_interp = np.interp(t_interp, t, x)# 时域信号绘图plt.figure(figsize=(10, 6))plt.subplot(2, 1, 1)# plt.plot(t_original, white_noise, 'b-', marker='o', label="Original Signal (fs=1000Hz)")plt.plot(t, x, 'b-', marker='o', linewidth=2,label="Original Signal (fs=1000Hz)")plt.plot(t_new, x_new, 'r--', marker='o', label="resampled Signal (fs=2000Hz)")plt.plot(t_interp, x_interp, 'y--', marker='o', label="interp Signal (fs=2000Hz)")plt.title("Original vs resampled Signal")plt.xlabel("Time [s]")plt.ylabel("Amplitude")plt.grid()plt.legend()# 频域信号绘图plt.subplot(2, 1, 2)frequencies, fft_amplitude = plot_frequency_response(x, fs)frequencies_new, fft_amplitude_new = plot_frequency_response(x_new, fs_new)frequencies_interp, fft_amplitude_interp = plot_frequency_response(x_interp, fs_new)plt.plot(frequencies, fft_amplitude, 'b-', marker='o', label="Original Signal (fs=1000Hz)")plt.plot(frequencies_new, fft_amplitude_new, 'r--', linewidth=2, marker='o', label="resampled Signal (fs=2000Hz)")plt.plot(frequencies_interp, fft_amplitude_interp, 'y--', linewidth=2, marker='o', label="interp Signal (fs=2000Hz)")# 细化刻度plt.gca().xaxis.set_major_locator(MultipleLocator(50))plt.xlabel('Frequency (Hz)')plt.ylabel('Amplitude')plt.grid(True)plt.title("Frequency Response")plt.legend()plt.tight_layout()plt.show()
感谢小伙伴的关注，后续将陆续为大家奉上数字数据处理原理及代码实现等内容....
如果您有场景、有算法...
何不考虑用我们的低成本采集器和免费软件进行二次开发创造自己的产品？
关注本号,联系我们，创造就开始了
千元机，8通道IEPE、ICP
可以快速部署算法及业务的免费软件平台
扫一扫
555
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
