实现傅里叶变换并作图的python代码

以下是Python中实现傅里叶变换并作图的代码示例：

import numpy as np
import matplotlib.pyplot as plt

# 生成需要处理的信号
t = np.linspace(0, 2 * np.pi, 1000)
y = 3 * np.sin(5 * t) + 2 * np.cos(10 * t) + np.random.rand(len(t))

# 计算信号的傅里叶变换
Y = np.fft.fft(y)

# 取前一半的幅值，排除直流分量和镜像部分
amplitude = 2 / len(t) * np.abs(Y[:len(t) // 2])

# 计算对应的频率
frequency = np.linspace(0, 1 / (2 * (t[1]-t[0])), len(t) // 2)

# 绘制图像
plt.figure()
plt.plot(frequency, amplitude)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude')
plt.show()

这段代码可以生成一个包含傅里叶变换结果的频谱图像，可以帮助我们可视化信号包含的频率信息。

相关问题

分别写一段根据桥梁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'])

傅里叶级数曲线拟合 Python

傅里叶级数是一种用正弦和余弦函数表示周期函数的方法。在Python中，可以使用scipy库中的fft函数来计算离散傅里叶变换（DFT），从而实现对周期信号的频域分析。可以利用傅里叶级数进行曲线拟合，具体步骤如下：

1. 采集数据并作图：首先，采集周期性数据，并将其用matplotlib库作图。

2. 计算傅里叶级数系数：通过离散傅里叶变换函数fft计算数据的傅里叶系数。

3. 选择拟合项：选择需要用于拟合的傅里叶级数项数，这取决于需要拟合的曲线复杂度。

4. 重构信号：使用所选的傅里叶级数项重构原始信号，并与原始信号进行比较。

5. 计算误差：根据重构信号和原始信号之间的误差，评估拟合的效果。

以下是一个简单的Python代码示例，用于实现傅里叶级数曲线拟合：

import numpy as np
import matplotlib.pyplot as plt

# 生成原始信号
t = np.linspace(-np.pi, np.pi, 201)
x = np.sin(t) + 0.5 * np.sin(3*t) + 0.2 * np.sin(5*t)

# 计算傅里叶系数
n = len(x)
k = np.arange(n)
T = n/2
frq = k/T
frq = frq[range(int(n/2))]
X = np.fft.fft(x)/n
X = X[range(int(n/2))]

# 选择拟合项
n_harm = 10
indexes = list(range(1, n_harm+1))
amplitudes = 2*np.abs(X[indexes])
phases = np.angle(X[indexes])

# 重构信号
t_recon = np.linspace(-np.pi, np.pi, 201)
x_recon = np.zeros_like(t_recon)
for i in range(n_harm):
    x_recon += amplitudes[i] * np.cos(i*t_recon + phases[i])

# 计算误差
error = np.mean((x - x_recon)**2)

# 绘制原始信号和重构信号
plt.plot(t, x, label='Original')
plt.plot(t_recon, x_recon, label='Reconstructed')
plt.legend()
plt.show()

# 输出相关问题