基于时频谱图像处理的海豚哨声信号提取算法

此方法针对复杂噪声环境中海豚哨声信号的提取问题，首先就是声学数据预处理，即将含有哨声信号的声学数据进行初步处理，如去除机械噪声和白噪声。然后进行时频谱(TFS)图像转换，使用短时傅里叶变换 (STFT)将声学数据转换为时频谱图像，这些图像能够展示哨声信号的时间-频率特性。以提高后续处理的准确性。通过对比信号和噪声的时频特性，可以更有效地进行哨声信号的提取。然后将复杂的时频谱图像转换为二值图像，利用形态学操作腐蚀噪声保留哨声。

一、哨声来源

原始的哨声信号来源：康奈尔大学麦考利图书馆收藏的中华白海豚哨声录制音频

二、主要环节

1.预处理

预处理采用了谱减法去除背景噪声。

2.生成时频谱

利用短时傅里叶变换（STFT）生成时频谱图像。

3.二值化

时频谱图像转换为二值图像。为了将 TFS 图像转换为二值图像，需要确定一个合适的阈值。

4.形态学操作

对二值图像进行形态学筛选时使用凸包查找Graham扫描算法及腐蚀操作。

5.逆变换

通过逆短时傅里叶变换（iSTFT）得到提取后时域表示信号。

三、Python代码实现

import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt
from scipy.signal import hann
import librosa  
# 假设声学数据存储在audio_data中
# 假设采样率为sample_rate
from scipy.io import wavfile
import cv2
# 读取音频文件
file_path = "your_audio_file.wav"
#audio_data, sample_rate = librosa.load("cwd.wav", sr=None)

# audio_data 是一个包含声音信号波形的一维numpy数组
# sample_rate 是这个声音信号的采样率，用来标识每秒采样的数据点数

def remove_noise(audio_data, sample_rate):
    # 使用短时傅里叶变换(STFT)将声学数据转换为时频谱图像
    _, _, Zxx = signal.stft(audio_data, fs=sample_rate, nperseg=256)

    # 计算频谱的幅度
    magnitude = np.abs(Zxx)

    # 计算过度减法因子
    alpha = 0.001

    # 计算幅度估计因子
    beta = 10

    # 应用谱减法去除噪声
    magnitude_subtracted = magnitude - alpha

    # 将幅度截断为非负值
    magnitude_subtracted = np.maximum(magnitude_subtracted, 0)

    # 进行幅度估计
    estimated_magnitude = magnitude_subtracted / (magnitude_subtracted + beta)

    # 将估计的幅度与原始相位相乘
    enhanced_magnitude = estimated_magnitude * np.exp(1j * np.angle(Zxx))

    # 通过逆短时傅里叶变换(iSTFT)增强声音信号
    _, enhanced_signal = signal.istft(enhanced_magnitude, fs=sample_rate)

    return enhanced_signal


# 步骤2: 使用短时傅里叶变换(STFT)将声学数据转换为时频谱图像
def stft(audio_data, sample_rate):
    f, t, Zxx = signal.stft(audio_data, fs=sample_rate, nperseg=256)
    return f, t, Zxx

# 步骤3: 将复杂的时频谱图像转换为二值图像，这里简单地选取一个阈值进行二值化
def binarize_spectrogram(Zxx):
    threshold = np.max(np.abs(Zxx)) *0.03 # 根据信号幅度选择阈值
    binary_image = np.abs(Zxx) > threshold
    return binary_image

# 步骤4: 利用图像处理算法将信号和噪声分离，这里简单地使用形态学操作进行分离
def separate_signal_and_noise(binary_image):
    from scipy.ndimage import binary_erosion, binary_dilation
    binary_image = binary_erosion(binary_image, iterations=1)
    #binary_image = binary_dilation(binary_image, iterations=1)
    return binary_image

def convexHull(binary_image):
    # 将二值图像数据类型转换为 uint8 类型
    binary_image_uint8 = np.uint8(binary_image * 255)

    # 使用 findContours 函数找到轮廓
    contours, _ = cv2.findContours(binary_image_uint8, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)

    # 创建一个与原始图像大小相同的空白图像
    mask = np.zeros_like(binary_image_uint8)

    # 绘制轮廓并找到凸包
    for contour in contours:
        hull = cv2.convexHull(contour)
        cv2.drawContours(mask, [hull], 0, (255, 255, 255), -1)

    # 将掩码图像转换为二值图像
    mask_binary = mask > 0
    result_image = binary_image_uint8 * mask_binary
    #cv2.imshow('Result Image', result_image)
    #cv2.waitKey(0)
    #cv2.destroyAllWindows()
    return result_image
def remove_vertical_noise(binary_image):
    # 定义腐蚀核
    kernel = np.ones((1, 3), np.uint8)  # 竖直核，可以根据需要调整大小

    # 应用腐蚀操作
    eroded_image = cv2.erode(binary_image, kernel, iterations=1)

    return eroded_image


# 步骤5: 通过逆短时傅里叶变换(iSTFT)增强哨声信号
def enhance_signal(Zxx, binary_image):
    enhanced_signal = signal.istft(Zxx * binary_image)[1]
    return enhanced_signal

# 绘制结果
def plot_results(audio_data, filtered_audio_data, f, t, Zxx, binary_image, enhanced_signal, sample_rate):
    plt.figure(figsize=(12, 10))
    plt.rcParams['font.sans-serif'] = ['SimHei']
    plt.subplot(3, 2, 1)
    plt.plot(audio_data)
    plt.title('原始音频')

    plt.subplot(3, 2, 2)
    plt.specgram(audio_data, Fs=sample_rate,)
    plt.title('原始时频谱')

    plt.subplot(3, 2, 3)
    plt.plot(filtered_audio_data)
    plt.title('预处理后音频')

    plt.subplot(3, 2, 4)
    plt.specgram(filtered_audio_data, Fs=sample_rate,)
    plt.title('预处理后时频谱')

    plt.subplot(3, 2, 5)
    plt.imshow(binary_image, aspect='auto', origin='lower',extent=[t.min(), t.max(), f.min(), f.max()])
    plt.title('二值图像')

    plt.subplot(3, 2, 6)
    plt.imshow(eroded_image, aspect='auto',origin='lower')
    plt.title('最终图像')
    
 

    plt.tight_layout()
    plt.show()



     #将增强后的音频保存为新文件
   

def load(path):
 audio_data, sample_rate = librosa.load(path, sr=None)
 return audio_data, sample_rate
    # 步骤1: 去除噪声
#path=input()
audio_data, sample_rate = load('cwd.wav')
filtered_audio_data = remove_noise(audio_data, sample_rate)

    # 步骤2: STFT
f, t, Zxx = stft(filtered_audio_data, sample_rate)

    # 步骤3: 二值化谱图
binary_image = binarize_spectrogram(Zxx)

    # 步骤4: 分离信号和噪声
binary_image = separate_signal_and_noise(binary_image)
result_image = convexHull(binary_image)
eroded_image = remove_vertical_noise(result_image)
    # 步骤5: 通过iSTFT增强信号
enhanced_signal = enhance_signal(Zxx, eroded_image)

    # 绘制结果
plot_results(audio_data, filtered_audio_data, f, t, Zxx, binary_image,enhanced_signal, sample_rate)
    #wavfile.write("enhanced_cwdaudio1.wav", sample_rate, enhanced_signal)

四、效果展示