np.fft.rfft处理混频AM波输出异常结果的问题排查

我是信号处理新手，尝试生成AM波后进行频率混频，并获取其FFT结果，代码如下：

import numpy as np
import matplotlib.pyplot as plt

# AM
fc = lambda t: np.sin(2 * np.pi * 10e6 * t)
fm = lambda t: 0.5 * np.sin(2 * np.pi * 3e3 * t)
f = lambda t: (1 + fm(t)) * fc(t)

# Frequency mix
fmix = lambda t: np.sin(2 * np.pi * 10e6 * t) # same to fc
f0 = lambda t: f(t) * fmix(t) + 1.5

N = 1024
FS = 2e4
T = (1.0 / FS) * N

t = np.linspace(0, T, N)
y = f0(t)
Y = np.fft.rfft(y)
Y_amp = abs(Y)
Y_amp_norm = Y_amp / N
freq = np.fft.rfftfreq(N, T / N)
idx = np.argsort(freq)

plt.plot(freq[idx], Y_amp_norm[idx])
plt.show()

我认为输出图像应分别在0Hz和3kHz处出现峰值，但实际图像的第二个峰值约在400Hz处，即使在FFT前添加低通滤波器，结果仍一致，请问整个流程哪里出错了？

核心错误点 #

1. 采样率严重不足导致混叠失真 #

你设置的载波频率是10MHz，但采样率FS只有20kHz，远远低于奈奎斯特采样定理要求的2倍载波频率（20MHz）。高频信号直接被折叠到低频段，这是峰值位置完全错误的根本原因——你看到的400Hz峰值其实是10MHz信号混叠后的虚假频率。

2. 混频后的信号处理逻辑冗余 #

根据三角恒等式，两个正弦信号相乘的结果是：
sin(A)*sin(B) = [cos(A-B) - cos(A+B)]/2
你的AM信号(1+fm(t))*sin(2π*10e6*t)和本振sin(2π*10e6*t)相乘后，应该得到：
0.5*(1+fm(t)) - 0.5*(1+fm(t))*cos(2π*20e6*t)
这里的0.5*(1+fm(t))就是包含3kHz分量的基带信号，而你额外加的+1.5属于冗余的直流偏移，不是核心问题，但20MHz的高频分量因采样率不足直接混叠，完全掩盖了真实的3kHz信号。

3. FFT幅度归一化方式有误 #

对于rfft（实信号FFT），幅度归一化应该区分处理：直流分量和Nyquist分量除以N，其余分量除以N/2，你直接除以N会导致幅度被缩小一半，不过这不是峰值位置错误的原因。

修正后的代码示例 #

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import butter, filtfilt

# 基础参数修正：满足奈奎斯特采样要求
fc = 10e6  # 载波频率10MHz
fm = 3e3   # 调制频率3kHz
FS = 25e6  # 采样率设为25MHz，远大于2*fc
N = 1024
dt = 1 / FS  # 采样间隔
t = np.arange(N) * dt  # 用采样间隔生成时间序列，避免linspace的精度问题

# 生成AM信号
carrier = np.sin(2 * np.pi * fc * t)
modulator = 0.5 * np.sin(2 * np.pi * fm * t)
am_signal = (1 + modulator) * carrier

# 混频操作
local_osc = np.sin(2 * np.pi * fc * t)
mixed_signal = am_signal * local_osc

# 低通滤波：滤除20MHz的高频分量，保留基带
b, a = butter(4, 5e3 / (FS/2), btype='low')  # 截止频率设为5kHz
filtered_signal = filtfilt(b, a, mixed_signal)

# FFT分析
Y = np.fft.rfft(filtered_signal)
Y_amp = np.abs(Y)
# rfft归一化：直流和Nyquist分量除以N，其余除以N/2
Y_amp_norm = np.zeros_like(Y_amp)
Y_amp_norm[0] = Y_amp[0] / N
Y_amp_norm[-1] = Y_amp[-1] / N
Y_amp_norm[1:-1] = Y_amp[1:-1] / (N/2)

freq = np.fft.rfftfreq(N, dt)

plt.plot(freq, Y_amp_norm)
plt.xlabel("频率(Hz)")
plt.ylabel("归一化幅度")
plt.xlim(0, 10e3)  # 聚焦0-10kHz频段，方便观察3kHz峰值
plt.grid(True)
plt.show()

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