简介面向西电通信工程学院学生及信号与系统课程学习者这份大作业资料围绕多普勒效应展开提供从原理分析、MATLAB 建模到实验验证的完整参考适合课程作业、专题复习与工程入门实践。压缩包共 6 个文件约 1.03MB包含 4 个 .m 脚本、1 份 PDF 实践指导与 1 份 DOCX 报告文档分别对应四道分题的求解、仿真代码和书面说明便于对照理解多普勒频移、频谱分析、信号采样与参数设置等知识点。资源已有 1005 人学习说明其在校内课程场景中具有一定参考价值。读者可据此梳理多普勒效应在雷达、声纳及无线通信中的表现参考脚本复现仿真流程、调参过程与结果验证并借助文档组织实验数据、分析结论与报告结构提升把理论知识转化为工程方案的能力。1. 从一段“听不出变化”的录音说起多普勒大作业到底在练什么手机贴着路边录一段电动车喇叭回放时几乎听不出音调变化把其中 0.2 秒的数据做 FFT主峰却从 2.0 kHz 挪到了 2.35 kHz。信号与系统课程大作业里的多普勒效应难的从来不是写公式而是把这段频移在离散采样、噪声和窗函数里稳定测出来。这套资源围绕多普勒效应大作业把回波建模、频移估计、时频分析、参数表和报告结论串成一条线适合通信、电子、自动化方向正在做课程设计的人也适合已经写完代码、但图表和结论对不上的读者。它解决的是“公式会推、代码跑不通、结果不敢写进报告”的尴尬不是单纯给一份仿真脚本。拿到手先别急着改参数把模型假设和验收指标看一遍后面每一步才知道为什么这么设。2. 多普勒频移公式落到代码回波建模、离散采样与频率分辨率先把模型钉死再谈代码。多普勒大作业最常见的翻车是把单程公式和双程公式混着用最后频移差了一倍图看起来还挺像。课堂演示常用声学信号接收端和发射端在同一侧时回波经历往返双程频移写成 (f_d2v f_0\cos\theta/c)如果只算声源移动、接收固定则是单程 (f_dv f_0\cos\theta/c)。电磁波场景把声速换成光速频移量会小很多参数表要跟着换不能直接套声学那组数。场景频移公式关键假设常见用途收发同端回波往返(f_d2v f_0\cos\theta/c)远场、匀速、单基地声呐、雷达演示单程声源移动(f_dv f_0\cos\theta/c)接收端静止课堂录音演示目标横向切过(\cos\theta\approx0)距离变化率接近零解释频移消失目标靠近(f_d0)正方向统一正频移目标远离(f_d0)正方向统一负频移我一般先固定一组能算出手算结果的参数(f_02000) Hz、(v30) m/s、(c340) m/s双程频移约 (2\times30\times2000/340352.9) Hz。这个数不大不小FFT 图上能看清又不至于超出低通带宽。采样率选 48 kHz不是越大越好而是给 (f_0f_d) 留出足够余量同时让后续降采样有空间。2.1 用复指数混频搭一个最小可跑仿真直接对实信号做 FFT峰值会落在 2.35 kHz 附近测频也能测但报告里不好解释“为什么射频附近还有一大堆分量”。更稳的做法是先把回波混频到基带把 352.9 Hz 的频移搬到零频附近再低通、加窗、测峰。下面这段代码是最小闭环跑完能打印理论频移和估计频移。import numpy as np from scipy.signal import butter, filtfilt fs 48000 # 采样率至少大于 2*(f0fd_max) T 0.5 # 观测时长决定频率分辨率 f0 2000 # 发射信号中心频率 v 30 # 目标径向速度 c 340 # 声速电磁波场景改成 3e8 fd_true 2 * v * f0 / c t np.arange(0, T, 1 / fs) tx np.cos(2 * np.pi * f0 * t) # 简化回波只保留频移相位和时延在作业里可另外补 rx np.cos(2 * np.pi * (f0 fd_true) * t) rng np.random.default_rng(7) rx rx 0.02 * rng.standard_normal(t.size) # 混频到基带把 f0 搬掉保留 fd lo np.exp(-1j * 2 * np.pi * f0 * t) base rx * lo # 低通滤掉 2*f0 附近的分量截止频率要高于最大频移 b, a butter(4, 500 / (fs / 2), btypelow) base_lp filtfilt(b, a, base) # 加 Hann 窗抑制泄漏再做 FFT N base_lp.size win np.hanning(N) spec np.fft.fftshift(np.fft.fft(base_lp * win)) freq np.fft.fftshift(np.fft.fftfreq(N, 1 / fs)) mask (freq 0) (freq 1000) fd_est freq[mask][np.argmax(np.abs(spec[mask]))] print(f理论频移: {fd_true:.2f} Hz) print(f估计频移: {fd_est:.2f} Hz) print(f频率分辨率: {fs / N:.2f} Hz)这段代码的逻辑是“混频—低通—加窗—FFT—取峰”。混频把 (f_0f_d) 搬到 (f_d)低通把 (2f_0f_d) 的高频分量压掉加窗是为了减少矩形窗带来的频谱泄漏。参数上fs48000给 2.35 kHz 留了足够余量T0.5对应频率分辨率 (1/T2) Hz测 352.9 Hz 够用低通截止 500 Hz 覆盖了当前频移但如果速度范围扩到 60 m/s频移约 705.9 Hz截止频率要跟着提到 900 Hz 以上否则峰值会被滤掉。filtfilt是零相位滤波避免群时延把峰值位置带偏如果作业要求实时处理改成lfilter并补偿时延。提示混频用的复本振频率必须和发射频率一致。复本振写成 (e^{-j2\pi f_0t}) 还是 (e^{j2\pi f_0t}) 只影响正负频率方向但报告里的“靠近为正”要统一否则频移符号会反。2.2 观测时长和采样率互相牵制很多人以为把 FFT 点数补到 65536频率分辨率就变细了结果峰值位置确实更“平滑”误差却没降。补零只是频域插值不增加信息真正决定能不能分开两个频率的是观测时长 (T)分辨率近似 (1/T)。下面的小脚本把同一套参数在不同 (T) 下跑一遍能直观看到误差怎么变。def measure_fd(fs, T, v30, f02000, c340, snr_db20): rng np.random.default_rng(0) t np.arange(0, T, 1 / fs) fd_true 2 * v * f0 / c rx np.cos(2 * np.pi * (f0 fd_true) * t) sigma 10 ** (-snr_db / 20) # 相对单位幅度信号的噪声标准差 rx rx sigma * rng.standard_normal(t.size) base rx * np.exp(-1j * 2 * np.pi * f0 * t) b, a butter(4, 500 / (fs / 2), btypelow) base_lp filtfilt(b, a, base) N base_lp.size spec np.fft.fftshift(np.fft.fft(base_lp * np.hanning(N))) freq np.fft.fftshift(np.fft.fftfreq(N, 1 / fs)) mask (freq 0) (freq 1000) return fd_true, freq[mask][np.argmax(np.abs(spec[mask]))] for T in [0.02, 0.05, 0.1, 0.2, 0.5, 1.0]: fd_true, fd_est measure_fd(48000, T) print(fT{T:4}s 理论{fd_true:7.2f}Hz 估计{fd_est:7.2f}Hz 误差{abs(fd_est-fd_true):6.2f}Hz)T0.02s 时分辨率只有 50 Hz352.9 Hz 的峰可能落到 350 Hz 或 400 Hz 的栅格上误差看着像算法不行其实是观测太短。T0.5s 分辨率 2 Hz误差通常能压到几赫兹以内。snr_db控制噪声强度sigma10**(-snr_db/20)是按单位幅度信号折算的如果发射信号幅度不是 1要先把噪声标准差按实际幅度缩放否则信噪比标称值和图上看到的不一致。2.3 窗函数和补零的边界加窗不是玄学 rectangular window 在频率不是 bin 整数倍时泄漏严重Hann 窗能把旁瓣压低但主瓣会变宽两个靠近的频率更不容易分开。多普勒大作业通常只有一个主峰Hann 窗是稳妥选择如果要做多目标分辨窗长和窗型要写进参数表不能只在代码里藏一个hann。补零可以用但要在报告里说明“补零只改善峰值插值不改变物理分辨率”否则验收时一问就露馅。还有一点容易漏实信号混频后是复数FFT 取正频率半边时不要重复计算能量。画频谱时我一般用20*np.log10(np.abs(spec)1e-12)加一个小量避免 log 零峰值搜索只在有效频段内做别把 0 Hz 直流和低通过渡带算进去。把这几步写清楚后面的时频图和速度曲线才有可信的起点。3. 时频图不是装饰用STFT把速度变化拆到每一帧单帧 FFT 给的是整段观测里的平均频移。目标如果先靠近再远离或者速度在变一个峰值会把两段运动糊在一起报告里写“平均速度 30 m/s”没问题但解释不了时频图上的斜线。短时傅里叶变换STFT就是把长信号切成小段每段做 FFT横轴时间、纵轴频率频移随时间怎么变一眼能看出来。它不提高频率分辨率只是把“什么时候有哪些频率”摊开所以窗长怎么选直接决定图能不能用。3.1 窗长、重叠和补零的三个参数STFT 的核心参数是nperseg、noverlap和nfft。nperseg是每段窗长对应时间分辨率约nperseg/fs频率分辨率约fs/npersegnoverlap决定帧移75% 重叠对应帧移nperseg/4图更连续但计算量更大nfft是 FFT 点数大于nperseg时做补零插值让脊线看起来更光滑但不增加真实分辨率。我一般先按“频率分辨率要小于最小待测频移的 1/3”估nperseg再检查时间分辨率能不能跟上速度变化。参数示例值物理含义调大后的影响nperseg1024窗内采样点数频率分辨率变好时间分辨率变差noverlap768相邻窗重叠点数帧移变小脊线更连续nfft4096FFT 点数频域插值更密真实分辨率不变windowhann窗型旁瓣低主瓣变宽3.2 一段能出脊线的STFT脚本下面这段代码接着前面的基带信号跑输出每一帧的峰值频率再换算成速度。画图部分只保留必要设置报告里可以把时频图和速度曲线上下对齐放。import numpy as np from scipy.signal import stft, medfilt fs 48000 f0 2000 c 340 # base_lp 来自上一节的混频低通结果 f, t_stft, Z stft( base_lp, fsfs, windowhann, nperseg1024, noverlap768, nfft4096, return_onesidedTrue, boundaryNone ) mag np.abs(Z) valid (f 100) (f 1000) # 避开直流和低通过渡带 idx np.argmax(mag[valid], axis0) f_peak f[valid][idx] # 频移换速度双程模型 v f_d*c/(2*f0) v_est f_peak * c / (2 * f0) v_smooth medfilt(v_est, kernel_size5) print(时间帧数:, t_stft.size) print(前 5 帧峰值频率:, np.round(f_peak[:5], 2)) print(前 5 帧速度估计:, np.round(v_smooth[:5], 2))nperseg1024在 48 kHz 下对应 21.3 ms 窗长频率分辨率约 46.9 Hznoverlap768是 75% 重叠帧移 5.33 msnfft4096把频域栅格插到 11.7 Hz脊线看起来更细但真正的分辨能力还是由 21.3 ms 窗长决定。valid掩掉 100 Hz 以下是因为低通过渡带和直流泄漏会把峰值拉到 0 Hz 附近尤其在小速度场景下更明显。medfilt的核大小 5 是经验值帧数少的时候别用太大否则速度突变会被抹平。3.3 从脊线反推速度曲线频移换速度时先把正负方向定死目标靠近、回波频率升高(f_d0)远离则 (f_d0)。如果 STFT 只取了正频率半边负频移会折到正半轴速度曲线会突然跳变。稳妥做法是保留复数 STFT在fftshift后的频率轴上取峰值或者对混频后的复基带做双边谱。报告里不要只放时频图最好把“时间—峰值频率—速度”三列数据导出成 CSV图和数据能对上验收时被追问也有据可查。时频图在报告里的作用不是好看而是解释“为什么单帧 FFT 的平均值不能代表全过程”。目标匀速时脊线是水平线靠近时在正频移区远离时在负频移区如果脊线断裂先查信噪比和重叠率再查峰值搜索有没有被直流或谐波带偏。把这三个判断写进图注比堆十张彩图更有说服力。4. 参数表不是拍脑袋从验收指标反推仿真配置多普勒大作业的验收指标通常不会写成“把图跑出来就行”而是隐含在问题里频移误差多大、速度误差多大、最小可测速度是多少、信噪比降到多少还能出脊线。参数表要从这些指标反推不是先写代码再补表。常见做法是先列一张指标表再算采样率、观测时长、低通截止和 FFT 点数最后用参数扫描验证边界。指标目标值反推参数备注速度范围060 m/s最大频移约 705.9 Hz声速 340 m/s、(f_02000) Hz频移误差10 Hz(T0.1) s分辨率 (1/T)速度误差1 m/s(T0.2) s(c/(2f_0)\approx0.085) m/s/Hz采样率大于 2 倍最高回波频率至少 5411 Hz实际取 48 kHz 留余量信噪比020 dB加噪脚本可扫报告里写清噪声定义按这张表(T0.2) s 时频率分辨率 5 Hz速度误差约 0.425 m/s满足 1 m/s 指标(T0.1) s 时分辨率 10 Hz速度误差约 0.85 m/s刚好卡在边缘。低通截止要高于最大频移 705.9 Hz取 900 Hz 或 1 kHz 更稳采样率 48 kHz 给 2.35 kHz 的回波留了足够余量也方便后续做降采样对比。4.1 用参数扫描找可用区间单点跑通不算数换几个速度和信噪比还能稳住才敢写进报告。下面这段扫描脚本把速度从 5 m/s 扫到 50 m/s、信噪比从 0 dB 扫到 20 dB输出平均绝对误差和最大绝对误差。import numpy as np from scipy.signal import butter, filtfilt def estimate_error(v, snr_db, fs48000, T0.2, f02000, c340): rng np.random.default_rng(42) t np.arange(0, T, 1 / fs) fd_true 2 * v * f0 / c rx np.cos(2 * np.pi * (f0 fd_true) * t) sigma 10 ** (-snr_db / 20) rx rx sigma * rng.standard_normal(t.size) base rx * np.exp(-1j * 2 * np.pi * f0 * t) b, a butter(4, 1000 / (fs / 2), btypelow) base_lp filtfilt(b, a, base) N base_lp.size spec np.fft.fftshift(np.fft.fft(base_lp * np.hanning(N))) freq np.fft.fftshift(np.fft.fftfreq(N, 1 / fs)) mask (freq 100) (freq 1000) fd_est freq[mask][np.argmax(np.abs(spec[mask]))] return abs(fd_est - fd_true) v_list np.array([5, 10, 20, 30, 40, 50]) for snr in [0, 5, 10, 20]: errs [estimate_error(v, snr) for v in v_list] print(fSNR{snr:2}dB MAE{np.mean(errs):6.2f}Hz MAX{np.max(errs):6.2f}Hz)snr_db按单位幅度信号折算噪声标准差和上一节保持一致T0.2是反推出来的观测时长低通截止放到 1000 Hz覆盖 50 m/s 对应的 588 Hz 频移。扫描结果里如果某个速度在 0 dB 下误差突然变大不要急着改算法先看峰值是不是落到了低通过渡带或者被噪声旁瓣顶掉。把这张表放进报告比只写“仿真效果良好”有说服力。4.2 加噪、去均值、加窗三个预处理噪声不是最后才加预处理也不是可有可无。回波里如果有直流偏置混频后 0 Hz 附近会出现一个大峰峰值搜索很容易被带偏先减去均值再做带通或低通能把直流压下去。加窗放在 FFT 之前不要放在混频之前否则会把有用信号边沿也削掉。import numpy as np from scipy.signal import butter, filtfilt # 假设 rx 是原始回波fs 已定义 x rx - np.mean(rx) # 去直流 # 带通下限压低频漂移上限压高频噪声 b, a butter(4, [100 / (fs / 2), 1500 / (fs / 2)], btypebandpass) x_bp filtfilt(b, a, x) # 加窗放在 FFT 前不要提前加 x_win x_bp * np.hanning(x_bp.size)带通下限 100 Hz 是为了避开直流和低频漂移但如果目标速度很低频移可能只有几十赫兹下限要跟着降否则会把真实频移滤掉。上限 1500 Hz 覆盖当前回波频率换到 60 m/s 场景要重新算。filtfilt零相位适合离线分析实时处理要换因果滤波器并补偿时延。加窗只影响频谱泄漏不改变信号能量报告里写清窗型和窗长即可。4.3 结果归档文件和参数要能对上参数表、脚本、图、结论四样东西必须能互相追溯。我一般按“场景_参数_版本”命名输出目录图里标注fs、T、nperseg、SNRCSV 里保留时间、峰值频率、速度三列。这样验收时被问到“这张图是哪个参数跑的”不用翻聊天记录。输出物命名示例必须包含时域图wave_v30_snr10.png采样率、时长、速度频谱图fft_v30_snr10.png窗型、FFT 点数、峰值频率时频图stft_v30_snr10.png窗长、重叠、nfft速度曲线speed_v30_snr10.csv时间、频率、速度参数表config_v30_snr10.json全部可调参数5. 排查清单多普勒大作业里最容易翻车的五个点这套流程跑通一次不难难的是换参数、换场景后还能复现。下面五条是我在类似大作业里见过最多的翻车现场每条按“现象—原因—解决”写照着排查能省不少来回。5.1 频域和混频上的三个坑现象峰值总在 0 Hz 附近速度估计接近零。原因通常是回波里带了直流偏置或者混频本振和发射频率不一致差出一个固定低频。先看时域均值是不是明显不为零再检查lo的频率是不是严格等于f0如果发射信号有起始相位混频后会有固定相偏但不会把峰值搬到 0 Hz。解决方法是去均值、加带通并把本振频率和发射频率写在同一个配置里别一个写 2000、一个写 2005。现象频移符号反了靠近的目标算出负速度。原因多半是正方向定义不统一公式里写了 (f_d2vf_0\cos\theta/c)代码里却用 (e^{j2\pi f_0t}) 混频正负频率对调。解决方法是先在报告里写死“靠近为正、远离为负”再用一个已知速度的仿真信号验证符号如果符号反了要么改混频指数符号要么在速度换算前乘 -1但只能改一处不能两处都改。现象补零到 65536 点后峰值更细误差却没降。原因是补零只做频域插值不增加观测信息真实分辨率还是 (1/T)。如果 (T0.02) s分辨率 50 Hz补零再多也只能在 50 Hz 栅格之间插值。解决方法是先加长观测时长或者提高信噪比让峰值更尖锐补零可以保留但报告里要写明“用于峰值插值不改变物理分辨率”别把补零当成提高精度的手段。5.2 时频图和交付上的两个坑现象STFT 脊线断成几截速度曲线跳变。原因可能是窗长太短、重叠不够或者峰值搜索被直流、谐波带偏。先看时频图在断裂处是不是信噪比很低再把noverlap提到 75% 以上峰值搜索加一个valid频段掩码避开 0 Hz 和低通过渡带。如果目标速度变化很快窗长要缩短但频率分辨率会变差这时要在报告里说明取舍不能既想要 1 Hz 分辨率又想要 1 ms 时间分辨率。现象报告里的参数和代码对不上图注写错采样率。原因是参数散落在脚本、笔记和图注里改了一处忘了另一处。解决方法是用一个config.json或params.py存全部可调参数脚本只读配置图标题从配置里拼字符串交稿前跑一遍“参数—图—结论”对照确认fs、T、nperseg、SNR四处一致。这个坑看着小验收时被追问一次就很被动。注意排查时先固定随机种子再改参数。噪声种子一变峰值位置可能差几个 bin容易把参数问题误判成算法问题。6. 验收前自检用一张对照表把参数、图和结论钉死最后一晚最容易出的问题不是代码跑不通而是图、参数、结论三张皮。我一般写一个很小的自检脚本读配置文件检查采样率是否满足奈奎斯特、观测时长是否够分辨率、低通截止是否覆盖最大频移再把关键数字打印出来和图注核对。import json import pathlib cfg json.loads(pathlib.Path(config.json).read_text()) fs cfg[fs] f0 cfg[f0] fd_max cfg[fd_max] T cfg[T] fc cfg[lowpass_hz] assert fs 2 * (f0 fd_max), 采样率不满足奈奎斯特 assert 1 / T cfg[freq_error_target_hz], 观测时长不够频率分辨率偏粗 assert fc fd_max * 1.1, 低通截止太接近最大频移 print(参数自检通过) print(f分辨率{1/T:.2f}Hz, 速度刻度{cfg[c]/(2*f0):.4f}m/s per Hz)这段脚本不复杂作用是把“我以为的参数”变成“文件里的参数”。fd_max按最大速度算好写进配置freq_error_target_hz是报告里承诺的频移误差lowpass_hz至少比最大频移高 10%。如果断言失败先改配置再跑图别手动改图注糊弄过去。对照项参数表图注结论检查方式采样率48 kHz48 kHz满足奈奎斯特自检脚本断言观测时长0.2 s0.2 s分辨率 5 Hz打印 (1/T)窗长10241024时间分辨率 21.3 ms图注核对信噪比10 dB10 dB峰值可测噪声脚本复跑速度误差1 m/s与表一致满足指标CSV 复算从那以后我每次交大作业前都会强制走一遍这张“参数—图—结论”对照表哪怕时间只剩半小时。先跑自检脚本再抽查两张图的图注最后把 CSV 里的峰值频率代回公式算一遍速度三处对得上才敢点提交。希望帮到你。本文还有配套的精品资源点击获取