分别写一段根据桥梁GNSS加速度数据去噪并获取桥梁频率,作图展示的matlab和python代码

Matlab代码：

%读取数据
data = load('GNSS_acc_data.txt');

%去除噪声
fs = 100; %采样率
fc = 10; %截止频率
[b,a] = butter(2,fc/(fs/2)); %设计2阶巴特沃斯滤波器
acc_filtered = filtfilt(b,a,data); %零相移滤波

%获取频率
N = length(acc_filtered); %数据长度
t = (0:N-1)/fs; %时间轴
acc_fft = abs(fft(acc_filtered)); %进行傅里叶变换
f = (0:N-1)*(fs/N); %频率轴
f_bridge = f(acc_fft == max(acc_fft)); %获取峰值频率

%作图展示
subplot(2,1,1);
plot(t,data);
xlabel('Time (s)');
ylabel('Acceleration (m/s^2)');
title('Original GNSS Acceleration Data');
subplot(2,1,2);
plot(t,acc_filtered);
xlabel('Time (s)');
ylabel('Acceleration (m/s^2)');
title('Filtered GNSS Acceleration Data');
figure;
plot(f,acc_fft);
xlabel('Frequency (Hz)');
ylabel('Amplitude');
title('FFT of GNSS Acceleration Data');
hold on;
plot(f_bridge,max(acc_fft),'ro');
legend('FFT','Bridge Frequency');

Python代码：

#读取数据
import numpy as np
data = np.loadtxt('GNSS_acc_data.txt')

#去除噪声
from scipy.signal import butter, filtfilt
fs = 100 #采样率
fc = 10 #截止频率
b, a = butter(2, fc/(fs/2), 'lowpass') #设计2阶巴特沃斯滤波器
acc_filtered = filtfilt(b, a, data) #零相移滤波

#获取频率
N = len(acc_filtered) #数据长度
t = np.arange(N)/fs #时间轴
acc_fft = np.abs(np.fft.fft(acc_filtered)) #进行傅里叶变换
f = np.arange(N)*(fs/N) #频率轴
f_bridge = f[acc_fft == max(acc_fft)] #获取峰值频率

#作图展示
import matplotlib.pyplot as plt
plt.subplot(2,1,1)
plt.plot(t,data)
plt.xlabel('Time (s)')
plt.ylabel('Acceleration (m/s^2)')
plt.title('Original GNSS Acceleration Data')
plt.subplot(2,1,2)
plt.plot(t,acc_filtered)
plt.xlabel('Time (s)')
plt.ylabel('Acceleration (m/s^2)')
plt.title('Filtered GNSS Acceleration Data')
plt.figure()
plt.plot(f,acc_fft)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude')
plt.title('FFT of GNSS Acceleration Data')
plt.plot(f_bridge,max(acc_fft),'ro')
plt.legend(['FFT','Bridge Frequency'])