laser_cal/frequency_filter.py

import pandas as pd
import matplotlib.pyplot as plt
import warnings
from pandas.errors import SettingWithCopyWarning
from scipy.signal import butter, filtfilt
from scipy.signal import welch
from scipy.fft import fft, fftfreq
import numpy as np


warnings.filterwarnings("ignore", category=SettingWithCopyWarning)
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False


def butter_lowpass_filter(data, cutoff, fs, order=5):
    nyq = 0.5 * fs  # Nyquist频率
    normal_cutoff = cutoff / nyq
    b, a = butter(order, normal_cutoff, btype='low', analog=False)
    y = filtfilt(b, a, data)
    return y


def butter_bandpass_filter(data, cutoff_low, cutoff_high, fs, order=5):
    nyq = 0.5 * fs  # Nyquist频率
    high_cutoff = cutoff_high / nyq
    low_cutoff = cutoff_low / nyq
    b, a = butter(order, [low_cutoff, high_cutoff], btype='band', analog=False)
    y = filtfilt(b, a, data)
    return y


def apply_fft(x, fs):
    n = len(x)
    t = 1.0 / fs
    fft_cof = fft(x - np.mean(x))  # Obtain FFT coefficients after removing the mean (DC component) from the signal.
    xf = fftfreq(n, t)[:n // 2]  # Taking positive spectrum only.
    # Multiply  abs(FFT coefficients) by 2 to compensate for positive spectrum and normalize by signal length.
    fft_positive = 2.0 / n * np.abs(fft_cof[0:n // 2])
    return xf, fft_positive


# Function to compute PSD
def compute_psd(signal, length, fs):
    seg_length = length / 10000  # Segment length.
    overlap = seg_length / 20  # overlap between segments (in number of samples)
    nfft_length = 2 ** 14  # FFT length
    frequencies, psd = welch(signal, fs=fs, window='han', nperseg=seg_length, noverlap=overlap, nfft=nfft_length)
    return frequencies, psd


def tower_cal(data, start_points, end_points, fs):

    """
    计算测量叶片数据的塔筒测量距离
    :param data: cycle_calculate计算完成后的数据。
    :param start_points: 所有每个周期开始点，叶片前缘突变点。
    :param end_points: 叶片后缘突变点。
    :param fs: 采样频率。
    """

    cutoff_low = 0.01 * fs  # 设置低通滤波截止频率
    combined_df_sorted = pd.concat([start_points, end_points]).sort_values(by='time')  # 合并DataFrame并按 'time' 列排序
    # 检查排序后的数据从end开始，start结束
    if combined_df_sorted.iloc[0].equals(start_points.iloc[0]):
        combined_df_sorted = combined_df_sorted.iloc[1:]
    if combined_df_sorted.iloc[-1].equals(end_points.iloc[-1]):
        combined_df_sorted = combined_df_sorted.iloc[:-1]
    combined_df_sorted.reset_index(drop=True, inplace=True)

    # 将 start_points 中的时间点转换为列表
    start_times = combined_df_sorted['time'].tolist()
    fft_positive_filters = []
    tower_dists = []
    # 遍历所有起始时间点
    for i in range(0, len(start_times), 2):

        # 获取当前起始和结束时间点
        start_time = start_times[i] + 0.01 * (start_times[i + 1] - start_times[i])
        end_time = start_times[i + 1] - 0.01 * (start_times[i + 1] - start_times[i])

        # 根据当前起始时间点和结束时间点对数据进行分段
        segment = data[(data['time'] > start_time) & (data['time'] <= end_time)]
        mean_point = segment['distance'].mean()
        tower_dists.append(mean_point)

    tower_dist = np.mean(tower_dists)

    return tower_dist


def process_fft(data, fs):

    """
    塔筒数据每组低通滤波，求平均得振动分析；
    叶片数据分组高通滤波
    :param data: tower_cal计算完成后的数据。
    :param fs: 采样频率。
    """

    cutoff_low = 0.001 * fs  # 设置低通滤波截止频率
    segment = data.head(int(len(data) * 0.95))

    # 确保segment的长度是2的n次幂
    desired_length = 2**((len(segment) - 1).bit_length() - 1)
    segment = segment.head(desired_length)
    segment.loc[:, 'distance_filtered'] = butter_lowpass_filter(segment['distance'], cutoff_low, fs)

    # 提取时间序列数据
    time_series_filter = segment['distance_filtered'].values  # 使用滤波后的距离数据
    xf_filter, fft_positive_filter = apply_fft(time_series_filter, fs)

    fft_positive_filters_truncated = fft_positive_filter[0:1000]
    xf_filter_truncated = xf_filter[0:1000]
    # 将 NumPy 数组转换为 Python 列表
    fft_y = fft_positive_filters_truncated.tolist()
    fft_y_scaled = [x * 1000 for x in fft_y]
    fft_y_scaled = [0] + fft_y_scaled
    fft_x = xf_filter_truncated.tolist()
    fft_x = [0] + fft_x
    fft_y_scaled = [0 if a_val < 0.1 else b_val for a_val, b_val in zip(fft_x, fft_y_scaled)]
    max_value = max(fft_y_scaled)
    max_index = fft_y_scaled.index(max_value)

    '''
    plt.plot(fft_x, fft_y_scaled, label='Filtered Signal')
    plt.xlabel('频率 (Hz)', fontsize=8)
    plt.ylabel('振幅 (m)', fontsize=8)
    plt.tick_params('x', labelsize=8)
    plt.tick_params('y', labelsize=8)
    plt.title('Frequency Spectrum', fontsize=12)
    plt.legend(fontsize=7)
    plt.grid(True)
    plt.savefig(f"filter_fft1.png", dpi=110, pil_kwargs={"icc_profile": False})
    plt.close()

    plt.plot(segment['time'], segment['distance_filtered'], label='Filtered Signal')
    plt.xlabel('时间 (s)', fontsize=8)
    plt.ylabel('振幅 (m)', fontsize=8)
    plt.tick_params('x', labelsize=8)
    plt.tick_params('y', labelsize=8)
    plt.title('Time-Domain Waveform', fontsize=12)
    plt.legend(fontsize=7)
    plt.grid(True)
    plt.savefig(f"filter_fft2.png", dpi=110, pil_kwargs={"icc_profile": False})
    plt.close()
    '''

    return segment, fft_x, fft_y_scaled, round(fft_x[max_index], 2), round(max_value, 2)