逆傅里叶变换分解合成正弦波时遇irfft参数错误求助

问题描述 #

修改FFT信号分解代码用于分解自制的3个正弦波叠加信号，尝试单独绘制各分量时，在irfft(ft_threshed)步骤触发如下错误：

raise ValueError(f"Invalid number of FFT data points ({n}) specified.")
ValueError: Invalid number of FFT data points (0) specified.

尝试处理ft_threshed为空的情况后，问题仍未解决，需排查原因。

相关代码 #

import numpy as np
from numpy.fft import fft, ifft, fftfreq, rfft, irfft, rfftfreq
from scipy import fftpack
import matplotlib.pyplot as plt
from matplotlib import style
import matplotlib.animation as anim
import time
style.use('dark_background')
fig,ax=plt.subplots(2,2,constrained_layout=False,figsize=(12,6))
t=[]
data=[]
start_time=time.time()
def decompose_into_n_signals(srs,n,t):
    fourier=rfft(srs)
    frequencies=rfftfreq(srs.size,d=2e-2/srs.size)
    out=[]
    for vals in np.array_split(frequencies,n):
        if vals.size==0:
            break
        ft_threshed=fourier.copy()
        ft_threshed[(vals.min()>frequencies)]=0
        ft_threshed[(vals.max()<frequencies)]=0
        # print(ft_threshed)
        if(ft_threshed.size==0):
            break
        out.append(irfft(ft_threshed))
    return out

def animate(self,t,data):
    t.append(time.time()-start_time)
    t=t[-10:]
    data=np.append(data,np.sin(2*np.array(t))+np.sin(3*np.array(t))+np.sin(4*np.array(t)))
    data=data[-10:]
    y=decompose_into_n_signals(data,3,t)

    ax[0][0].clear()
    ax[0][0].set_title('Unfiltered Data')
    ax[0][0].set_ylim(-3,3)
    ax[0][0].set_ylabel('Signal')
    ax[0][0].set_xlabel('Time')
    ax[0][0].plot(t,data)

    ax[0][1].clear()
    ax[0][1].set_title('Signal 1')
    ax[0][1].set_ylim(-3,3)
    ax[0][1].set_ylabel('Signal')
    ax[0][1].set_xlabel('Time')
    # print(len(t),(y[0]).size)
    # print(y[0])
    ax[0][1].plot(t,y[0])

    ax[1][0].clear()
    ax[1][0].set_title('Signal 2')
    ax[1][0].set_ylim(-3,3)
    ax[1][0].set_ylabel('Signal')
    ax[1][0].set_xlabel('Time')
    ax[1][0].plot(t,y[1])

    ax[1][1].clear()
    ax[1][1].set_title('Signal 3')
    ax[1][1].set_ylim(-3,3)
    ax[1][1].set_ylabel('Signal')
    ax[1][1].set_xlabel('Time')
    ax[1][1].plot(t,y[2])

ani=anim.FuncAnimation(fig,animate,fargs=(t,data),interval=10,frames=100)
fig.tight_layout(pad=5.0)
plt.show()

问题原因 #

1. 频率分割逻辑错误：np.array_split(frequencies, 3)是对频率数组做均匀拆分，但目标信号的频率是2、3、4，这种拆分方式无法保证每个频率被分到单独区间，导致部分区间内没有有效FFT分量，ft_threshed被置0后无有效数据点，触发irfft报错。
2. 采样间隔计算错误：rfftfreq的d参数是采样间隔，代码中2e-2/srs.size完全不符合实际采样间隔（实际应为相邻时间点的差值），导致频率轴计算错误，无法正确匹配信号频率。
3. 初始数据量不足：动画迭代初期t和data长度极短（如第一次迭代仅1个数据点），rfft生成的频率点极少，分割后更容易出现空分量。

修复方案 #

修改要点 #

• 修正rfftfreq的d参数，使用实际采样间隔
• 替换频率分割逻辑，改为根据目标频率（2、3、4）精准筛选FFT分量
• 初始化时设置足够的初始数据长度，避免迭代初期数据量不足

修改后的完整代码 #

import numpy as np
from numpy.fft import rfft, irfft, rfftfreq
import matplotlib.pyplot as plt
from matplotlib import style
import matplotlib.animation as anim
import time

style.use('dark_background')
fig, ax = plt.subplots(2, 2, constrained_layout=False, figsize=(12, 6))

# 初始化足够长度的时间和数据，避免初期数据量不足
t = np.linspace(0, 0.2, 10).tolist()
data = (np.sin(2 * np.array(t)) + np.sin(3 * np.array(t)) + np.sin(4 * np.array(t))).tolist()
start_time = time.time()

def decompose_into_target_signals(srs, target_freqs, t):
    fourier = rfft(srs)
    # 计算正确的采样间隔：相邻时间点的平均差值
    dt = np.diff(t).mean()
    frequencies = rfftfreq(srs.size, d=dt)
    out = []
    
    for freq in target_freqs:
        ft_threshed = fourier.copy()
        # 给目标频率设置一个小的容差，匹配FFT的频率分辨率
        tol = 0.1
        mask = (np.abs(frequencies - freq) > tol)
        ft_threshed[mask] = 0
        # 检查是否有非零分量
        if np.all(ft_threshed == 0):
            # 如果没有有效分量，生成和输入同长度的零数组
            out.append(np.zeros_like(srs))
        else:
            out.append(irfft(ft_threshed))
    return out

def animate(self, t, data):
    # 添加新的时间点，保持固定采样间隔（0.02秒）
    new_t = time.time() - start_time
    t.append(new_t)
    t = t[-10:]
    
    # 生成新的信号数据
    new_data = np.sin(2 * new_t) + np.sin(3 * new_t) + np.sin(4 * new_t)
    data.append(new_data)
    data = data[-10:]
    
    # 分解为目标频率的分量（2、3、4）
    y = decompose_into_target_signals(np.array(data), [2, 3, 4], t)

    # 绘制原始信号
    ax[0][0].clear()
    ax[0][0].set_title('Unfiltered Data')
    ax[0][0].set_ylim(-3, 3)
    ax[0][0].set_ylabel('Signal')
    ax[0][0].set_xlabel('Time')
    ax[0][0].plot(t, data)

    # 绘制分量1（频率2）
    ax[0][1].clear()
    ax[0][1].set_title('Signal 1 (Freq=2)')
    ax[0][1].set_ylim(-1.5, 1.5)
    ax[0][1].set_ylabel('Signal')
    ax[0][1].set_xlabel('Time')
    ax[0][1].plot(t, y[0])

    # 绘制分量2（频率3）
    ax[1][0].clear()
    ax[1][0].set_title('Signal 2 (Freq=3)')
    ax[1][0].set_ylim(-1.5, 1.5)
    ax[1][0].set_ylabel('Signal')
    ax[1][0].set_xlabel('Time')
    ax[1][0].plot(t, y[1])

    # 绘制分量3（频率4）
    ax[1][1].clear()
    ax[1][1].set_title('Signal 3 (Freq=4)')
    ax[1][1].set_ylim(-1.5, 1.5)
    ax[1][1].set_ylabel('Signal')
    ax[1][1].set_xlabel('Time')
    ax[1][1].plot(t, y[2])

ani = anim.FuncAnimation(fig, animate, fargs=(t, data), interval=20, frames=100)
fig.tight_layout(pad=5.0)
plt.show()

内容的提问来源于stack exchange，提问作者Ritwik Das