Python设计FIR滤波器和正交移相器

YFC0001

已于 2023-12-25 20:03:44 修改

阅读量715

点赞数 6

文章标签： python numpy 开发语言

于 2023-12-23 23:15:20 首次发布

本文链接：https://blog.csdn.net/m0_67424508/article/details/135176338

版权

作者所见题目：

1.设计根据希尔伯特变换的定义，在一个时域上写一个长度为N的希尔伯特变换器的系数，然后利用numpy库生成N阶的矩形窗，海宁窗，布莱克曼窗，汉明窗。两者相乘得到冲激响应。对于输入信号，首先时延N/2，得到正交移相器的实部，其次信号过希尔伯特系统得到正交移相器的虚部，返回复信号正交移相器。

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import matplotlib
matplotlib.use("Qt5Agg")  # 修改配置的后端 backend
# 汉字字体，优先使用楷体，找不到则使用黑体
plt.rcParams['font.sans-serif'] = ['Kaitt', 'SimHei']
# 正常显示负号
plt.rcParams['axes.unicode_minus'] = False
import matplotlib.pyplot as plt
import Myhilbert
import matplotlib
matplotlib.use("Qt5Agg")  # 修改配置的后端 backend
import soundfile as sf
Win=['hann','hamming','blackman','rect']
filename1=['music1u-32.wav','music2u-44.wav','music3e-44.wav','music4e-32.wav','music5u-32.wav']
i=1
for fn in filename1:
    # 读取音频文件
    audio, samplerate = sf.read(fn)
    for win in Win:
        orthogonality_data, h = Myhilbert.hilbert(311, win, audio)
        # 图绘制原信号与滤波后信号的时域
        plt.figure(i)
        plt.title(f'窗为：{win},音频为：{fn}')
        plt.plot(np.arange(len(audio)), audio, label=f'原信号波形')
        plt.plot(np.arange(len(orthogonality_data.imag)), orthogonality_data.imag,label=f'信号过正交滤波器')
        # 添加图例
        plt.legend()
        i=i+1

以下是音频信号时域与移相.py文件的代码。其主体是两个for循环，实现的功能是把同一个音频过由'hann','hamming','blackman','rect'所设计的正交移相器的时域波形连同其原时域波形画在同一张图中。

import numpy as np
def hilbert(N,win, data):#N是正交移相器的阶数,win-窗类型，data-音频
    Nw=N
    # 窗函数
    if( win=='rect'):
        window=win_rect = np.ones(Nw)  # 矩形窗函数
    if win=='hann':
        window= np.hanning(Nw)  # 海宁窗函数
    if win=='blackman':
        window = np.blackman(Nw)  # 布莱克曼窗函数
    if win=='hamming':
        window = np.hamming(Nw)  # 汉明窗函数
    # 正交移相器
    wins = np.concatenate((window, np.zeros(N - Nw)))
    h = np.zeros(N)
    m = int((N - 1) / 2)
    h[0:N:2] = 2 / (((np.arange(0, N, 2)) - m) * np.pi)
    h = h * wins
    signal_delayed = np.roll(data, N // 2)
    i = np.convolve(signal_delayed, h, mode='same')
    result = signal_delayed + 1j * i
    return result ,h

3．以下是分析性能.py文件的代码。其主体是hilbert函数和一个循环体。hilbert函数返回时域参数，循环体将同一窗的不同阶数正交移相器的的频域增益波形画出。

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import matplotlib
matplotlib.use("Qt5Agg")  # 修改配置的后端 backend
from scipy import signal
# import Myhilbert
import soundfile as sf
matplotlib.use("Qt5Agg")  # 修改配置的后端 backend
# 汉字字体，优先使用楷体，找不到则使用黑体
plt.rcParams['font.sans-serif'] = ['Kaitt', 'SimHei']
# 正常显示负号
import numpy as np
def hilbert(N,win, data):#N是正交移相器的阶数,win-窗类型，data-音频
    Nw=N
    # 窗函数
    if( win=='rect'):
        window=win_rect = np.ones(Nw)  # 矩形窗函数
    if win=='hann':
        window= np.hanning(Nw)  # 海宁窗函数
    if win=='blackman':
        window = np.blackman(Nw)  # 布莱克曼窗函数
    if win=='hamming':
        window = np.hamming(Nw)  # 汉明窗函数
    # 正交移相器
    wins = np.concatenate((window, np.zeros(N - Nw)))
    h = np.zeros(N)
    m = int((N - 1) / 2)
    h[0:N:2] = 2 / (((np.arange(0, N, 2)) - m) * np.pi)
    h = h * wins
    signal_delayed = np.roll(data, N // 2)
    i = np.convolve(signal_delayed, h, mode='same')
    result = signal_delayed + 1j * i
    return h
audio, samplerate = sf.read('music1u-32.wav')
plt.rcParams['axes.unicode_minus'] = False
Win=['hann','hamming','blackman','rect']
n=[51,100,151,200]
i=1
for win in Win:
    j = 1
    for N in n:
        # 绘图
        # 绘制图
        h = hilbert(N, win, audio)
        plt.figure(i)
        plt.subplot(2,2,j)
        h=hilbert(N,win,audio)
        h=np.abs(np.fft.fft(h))
        h=np.array(20*np.log10(h))
        plt.plot(np.arange(len(h)),h)
        plt.title(f'{win}窗正相滤波器，阶数：{N}')
        plt.ylabel('增益(dB)')
        j=j+1
    i=i+1
plt.show()

4.以下是滤波生成音频.py文件的代码。其主体是三个for循环，其实现的作用是创建一个f"D:\Thr_1\FIR_github\wav_{fn}_{j}"文件，用来存放对应的音乐。以阶数，窗为元素命名对应的音乐。J表示迭代次数，首先在阶数为51，100，151，200时进行滤波，对音乐的音质进行分析，再进行第二次乃至第三次第四次滤波。

import soundfile as sf
import os
import scipy.io.wavfile as wav
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import matplotlib
from scipy import signal
import Myhilbert
matplotlib.use("Qt5Agg")  # 修改配置的后端 backend
Win = ['hann', 'hamming', 'blackman', 'rect']
n = [51, 100, 151, 200]



# 指定文件夹路径
wave_name = ['music1u-32.wav',
             'music2u-44.wav',
             'music3e-44.wav',
             'music4e-32.wav',
             'music5u-32.wav']
j = 1  # 迭代次数
# n = [11, 21, 31, 41]#第一次迭代阶数
n=[51,100,151,200]#第二次迭代次数

for fn in wave_name:
    folder_path = f"D:\Thr_1\FIR_github\wav_{fn}_{j}"
    # 使用os模块中的mkdir()函数创建文件夹
    os.mkdir(folder_path)
    for N in n:
        for win in Win:
            # 读取音频文件
            samplerate, audio = wav.read(fn)
            orthogonality_data, h = Myhilbert.hilbert(311, win, audio)
            # 设置滤波器参数
            cutoff_freq = 7000  # 设置截止频率为7kHz
            nyquist_freq = 0.5 * samplerate
            num_taps = N  # FIR滤波器阶数
            taps = signal.firwin(num_taps, cutoff_freq /
                                 nyquist_freq, window=win)  # 创建窗的FIR滤波器系数
            # 对音频数据进行滤波
            filtered_audio = signal.convolve(
                orthogonality_data.imag, taps, mode='same')
            # 将音频数据转换为整数类型
            filtered_audio = filtered_audio.astype(np.int16)
            # 保存文件地址
            output_file = f"D:\Thr_1\FIR_github\wav_{fn}_{j}\{win}_阶数{N}_fir.wav"
            # 保存滤波后的音频文件
            wav.write(output_file, samplerate, filtered_audio)

5.下面是文件绘制语谱图代码.py的代码：这段代码主要实现了以下功能：

导入必要的库：numpy、matplotlib、pandas、Myhilbert、librosa、scipy.signal。

设置窗函数的选择和阶数列表。

读取音频文件并转换为幅度谱。

使用自定义的Myhilbert.hilbert()函数进行Hilbert变换，得到正交数据和滤波器系数。

设置滤波器的参数，并使用signal.firwin()函数生成FIR滤波器系数。

对正交数据进行滤波，得到滤波后的数据。

将滤波后的数据进行傅里叶变换，并转换为幅度谱。

绘制原信号和滤波后信号的语谱图。

整体来说，这段代码实现了对音频信号进行滤波，并绘制出原信号和滤波后信号的语谱图。

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import matplotlib
import Myhilbert
matplotlib.use("Qt5Agg")  # 修改配置的后端 backend
import librosa
import librosa.display
from scipy import signal
Win = ['hann', 'hamming', 'blackman', 'rect']
win=Win[0]
n = [11, 21, 31, 41,51,61,100,151,200]#第一次迭代阶数
# n = [11, 21, 31, 41]#第一次迭代阶数
# n=[51,100,151,200]#第二次迭代次数
N=n[3]
# 读取音频文件并转换为幅度谱
audio_file = 'music1u-32.wav'
y, sr = librosa.load(audio_file,sr=32000)
D = librosa.amplitude_to_db(np.abs(librosa.stft(y)), ref=np.max)
samplerate, audio=sr,y
audio = np.frombuffer(audio, dtype=np.int16)  # 将数据转换为numpy数组
orthogonality_data, h = Myhilbert.hilbert(311, win, audio)


# 设置滤波器参数
cutoff_freq = 7000  # 设置截止频率为7kHz
nyquist_freq = 0.5 * samplerate
num_taps = N  # FIR滤波器阶数
taps = signal.firwin(num_taps, cutoff_freq /
                     nyquist_freq, window=win)  # 创建窗的FIR滤波器系数
# 对音频数据进行滤波
filtered_audio = signal.convolve(
        orthogonality_data.imag, taps, mode='same')
fft_filtered_audio=librosa.amplitude_to_db(np.abs(librosa.stft(filtered_audio)), ref=np.max)
# fft_filtered_audio = np.squeeze(fft_filtered_audio)
#
# plt.figure()
# plt.plot(fft_filtered_audio.shape)
# plt.show()
# 绘制语谱图
plt.figure(figsize=(10, 8))
plt.subplot(2,1,1)
librosa.display.specshow(D, sr=sr, x_axis='time', y_axis='linear')
plt.colorbar(format='%+2.0f dB')
plt.title('原信号语谱图')
# plt.tight_layout()

plt.subplot(2,1,2)
librosa.display.specshow(fft_filtered_audio, sr=sr, x_axis='time', y_axis='linear')
plt.colorbar(format='%+2.0f dB')
plt.title(f'滤波后的语谱图,{N}阶,{win}窗')
plt.tight_layout()

# 显示语谱图
plt.show()
#
# # 将语谱图保存到文件
# # plt.savefig('spectrogram.png')

6.以下是研究不同采样频率对滤波器阶数的影响.py的代码：
这段代码是用于对音频信号进行滤波和频谱分析的。它使用了numpy、matplotlib、pandas和scipy等库。首先，通过导入相应的库，包括numpy、matplotlib.pyplot、pandas、matplotlib、wave和scipy.io.wavfile等。
然后，定义了一个窗口函数的列表Win，用于选择不同的窗口函数。接着，指定了要处理的音频文件路径wave_name和阶数n。
接下来，使用wave库打开音频文件，并读取音频数据。然后，将数据转换为numpy数组。接下来，对音频数据进行快速傅里叶变换（FFT），得到频谱数据。然后，计算频率轴，并进行频谱平移，使频谱图以0为中心频率。
然后，调用自定义的Myhilbert模块，使用正交希尔伯特变换对音频信号进行处理，得到正交信号和正交移相器系数。通过卷积操作，对音频信号进行滤波，得到滤波后的信号。接下来，使用scipy.signal库设置滤波器参数，并根据参数创建FIR滤波器系数。
然后，使用卷积操作对正交信号进行滤波，得到最终滤波后的信号。
最后，使用matplotlib.pyplot库绘制滤波后的信号的频谱图。整个过程会遍历不同的音频文件和阶数，每次都会绘制相应的频谱图。最后调用plt.show()来显示所有的频谱图。

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import matplotlib
import wave
import scipy.io.wavfile as wav

import Myhilbert

matplotlib.use("Qt5Agg")  # 修改配置的后端 backend
Win = ['hann', 'hamming', 'blackman', 'rect']
# 指定文件夹路径
wave_name = ['music1u-32.wav',
             'music2u-44.wav',
             'music3e-44.wav',
             'music4e-32.wav',
             'music5u-32.wav']
# wave_name=['music2u-44.wav']
wave_name = ['music1u-32.wav',
             'music2u-44.wav']
# n = [11, 21, 31, 41,51,61,100,151,200]#第一次迭代阶数
n = [11, 21, 31, 41]#第一次迭代阶数
# n=[51,100,151,200]#第二次迭代次数
n=[31]

import wave
# fn1=wave_name[0]
win=Win[1]

for fn1 in wave_name:
    for N in n:
        with wave.open(fn1, 'r') as wf:
            # 读取音频数据
            samplerate, audio = wav.read(fn1)
            audio = np.frombuffer(audio, dtype=np.int16)  # 将数据转换为numpy数组
        fft_audio=np.abs(np.fft.fftshift(np.fft.fft(audio)/(N/2)))
        M=len(fft_audio)
        f=2*np.pi*samplerate*(np.arange(M))/M
        f=f-max(f)/2
        orthogonality_data, h = Myhilbert.hilbert(311, win, audio)
        filtered_audio = np.convolve(h,audio,mode='same')
        fft_filtered_audio=np.abs(np.fft.fftshift(np.fft.fft(filtered_audio)/(N/2)))
        M=len(fft_filtered_audio)
        f=2*np.pi*samplerate*(np.arange(M))/M
        f=f-max(f)/2
       
        from scipy import signal
        # 设置滤波器参数
        cutoff_freq = 7000  # 设置截止频率为7kHz
        nyquist_freq = 0.5 * samplerate
        num_taps = N  # FIR滤波器阶数
        taps = signal.firwin(num_taps, cutoff_freq /
                             nyquist_freq, window=win)  # 创建窗的FIR滤波器系数
        # 对音频数据进行滤波
        filtered_audio = signal.convolve(
            orthogonality_data.imag, taps, mode='same')
        fft_filtered_audio=np.abs(np.fft.fftshift(np.fft.fft(filtered_audio)/(N/2)))
        M=len(fft_filtered_audio)
        f=2*np.pi*samplerate*(np.arange(M))/M
        f=f-max(f)/2
        plt.figure()
        plt.plot(f,fft_filtered_audio)
        plt.title(f'{fn1}信号过{win},{N}阶滤波器后的幅频')
plt.show()

有问题的欢迎评论留言。