公司动态
2024数维杯国际赛A题题解
数模没有标准答案仅供参考。模型的建立与求解题目预处理根据附件中的信息每个问题需要使用的表格数据不同我们首先将表格数据分成由问题一~问题四的四个表格信息方便于后续代码编写的计算。我们先对表达式进行变换转化为计算噪声特征z(t)的表达式问题一将频率、振幅和相位的值分别带入表达式频率f030´106Hz振幅A4相位j45°。问题二将振幅和相位的值分别带入表达式振幅A2相位j0°。将已知量带入表达式便于后续问题一、二的计算。问题1的模型建立与求解我们根据提供的表格信息计算噪声特征通过图表的方式直观感受对飞行周期1中接受数据的噪声特征进行分析1 绘制噪声特征值图像1根据预处理对问题一公式的处理计算得到的噪声特征值z_1(t)分析噪声的特征可以先通过自相关函数工具绘制出归一化后噪声信号的自相关函数描述信号在不同时间延迟下的相似度推测是否存在异常值。自相关函数的定义对于一个给定的时间序列xt其自相关函数rk定义为xt其中k—— 时间间隔滞后—— 时间序列的平均值通过绘制噪声的自相关变化图像我们可以发现信号或时间序列与自身在不同时间间隔时延下的相似程度如图1。图1自相关图像根据自相关图像分析在零延迟处的自相关值通常代表了信号的最大自相关归一化后这个值固定为1其他延迟位置的自相关值则反映了信号在这些延迟下相对最大自相关性的相对强度。2.通过绘制噪声振幅随时间变化图像如图2总结并分析噪声的变化特征。图2噪声随时间振幅的变化2 分析噪声特征值变化通过噪声值的自相关图像变化以及振幅随时间变化图像。总结到噪声特征值变化具有周期性噪声变化周期T约为0.0002s。根据实际情况噪声的特征值反映的就是收到信号的变化周期。根据噪声判断噪声的类型周期1的情况应该是色噪声。问题2的模型建立与求解1 问题模型的建立我们首先通过绘制图像观察噪声信号时域变化情况如图3图3周期2噪声信号时域变化情况如图3表示噪声信号的变化是平坦的不具有明显的峰值变化所以需要对信号进行变化处理。已知预处理得到的噪声特征值z_2(t)解决如何得到信号无噪声部分的特征值通过傅里叶变换将信号的特征时域转换到频域。通过描绘频谱图观察噪声部分的频谱变化。变换公式表示X[k]其中x[N]——时间序列信号N——数据点的数量X[k]——频率域上的值i——虚数单位在通过得到的频谱图找到振幅变化的最大点依据最大点计算出频谱中所对应的最大频率。其中k——从0到N-1的整数索引N——是样本点数Ts——采样时间间隔2飞行周期2频率估计计算根据上述问题模型的建立首先将通过傅里叶变换将接收到的信号从时域转换到频域绘制噪声信号的频谱图如图4。图4噪声信号的频谱图通过观察频谱图线性幅度变化找到了频谱中幅度最大的点。依据最大频率计算公式得到最大的频率的估计值3结果问题二依据噪声信号在时域变化中特征不明显我们将时域依据傅里叶公式转化到频域层面如图4呈现出来的峰值点依据峰值点估计出噪声信号的频率值。频率值即fmax3.82×106Hz。问题3的模型建立与求解1问题模型的建立结合实际情况振幅和相位同样是未知的结合前面的问题思路问题3主要解决的是噪声信号的振幅和相位的值。通过估计出的值再进行对飞行周期3情况下噪声信号频率值的估计。如何才能解决噪声信号振幅和相位值的问题众说周知噪声信号存在一种白噪声的情况这时它的均值为零。也就不用考虑振幅和相位的值。白噪声即具有零均值和方差已知。那如何解决频率问题根据时域的图像特征如图5与问题2相同我们使用问题二的方法计算通过傅里叶变换再通过使用最大似然估计频率f_3估计问题。图5周期3噪声信号时域变化情况例如问题2同样是找到振幅最大值点的频率。对频谱幅度取对数以更好地展示动态范围大的数据其中V是一个常数防止出现log(0)的情况。接下来找到峰值点对应的频率找到平滑后振幅的最大值及其频率2飞行周期3频率估计计算根据模型建立我们使用的是问题2的模型思路寻找振幅的最大值点计算出目标频率。通过我们绘制出结果图像如图6直观展示频率值。图6噪声信号的频谱图及结果4结果根据分析问题结合实际情况分析噪声信号的振幅和相位的值未知。直接计算需要考虑的因素居多我们引出噪声信号中白噪声的情况。在问题2模型的前提下我们加入最大似然估计对噪声信号进行估计计算。计算结果即3.5×107Hz。问题4的模型建立与求解1间歇接收条件下模型建立为了解决多信号情况下信号之间的干扰我们提出了间歇接收的方法进行信号采集。因为信号存在间歇性所以信号信息量将对较少因此频率估计变得困难。首先理解间歇接收模式间歇接收模式是信号在不连续的时间间隔内进行接收。是对某一时刻信号信息的采集。假设信号的接收过程中存在间歇性即存在某一时间点的信号样本。已知信号接收计算表达式x(t)Asin(2pf0tj)z(t)其中z(t)——噪声部分属于是高斯白噪声。时间不均匀进行采集我们需要对时间间隔进行计算即由于不连续采集信号信息较少。我们需要解决这个问题通过傅里叶变换补充时间不均匀情况下采集的信号信息。可以通过均值对不均匀处进行补充。通过对信号采样进行拟合最小化减少误差。因为会遇到采集信号变化较大可以实用最小二乘法对采样点的权重进行处理。权重的大小取决于采样间隔的大小。综上所述我们需要修正傅里叶变化的模型建立。在问题2的模型前提下我们加入不均匀时间间隔下被采样的情况。引入时间间隔。使用插值对信号情况进行补充利用信号x(α),其中ɑ是等时间间隔的时间点。2根据分析模型求解我们通过修正的傅里叶变换模型根据时域的变化依据不同的时间间隔。情况一我们删除很小的时间间隔留下时间间隔适中的信号情况计算出平均的时间间隔时间间隔情况如图所示。如图7时间间隔情况一。图7情况一平均时间间隔3.17×10-7S间隔较小我们舍去该情况。情况二我们保留较小的时间间隔通过计算计算出平均时间间隔时间间隔情况如图所示。如图8时间间隔情况二。如图8情况二平均时间间隔7.03×10-7S间隔较大我们保留该情况。我们使用情况二对问题的频率进行求解。求解结果如图9。直观看到最大值点以及最大频率的情况。图9 频率结果图3结果为了解决问题四提出的间歇性的实际情况我们通过分析时间序列的特征对信号特征的选取对信号特征的计算筛选出合适的平均时间间隔再利用修正傅里叶变化的模型对问题进行求解。通过计算振幅最大点进一步对频率结果进行估计估计出的目标结果频率即159694.319 Hz。代码问题一import numpy as np import pandas as pd import matplotlib.pyplot as plt from scipy.fft import fft, fftfreq from scipy.signal import correlate # 1. 读取 Excel 文件中的数据 file_path rD:\A题\pythonProject1\第一问.xlsx # 假设数据存储在第一个工作表中第一列是时间第二列是接收到的信号值 data pd.read_excel(file_path, sheet_name0) # 提取时间和接收到的信号值 time data[Received Signal Time].values x_t data[Received Signal Value].values # 2. 已知信号部分 (频率f030MHz幅度A4阶段phi45度) A 4 f0 30e6 # 30 MHz phi np.pi / 4 # 45 degrees # 重建已知信号部分 x_signal A * np.sin(2 * np.pi * f0 * time phi) # 3. 计算噪声部分 z(t) x(t) - x_signal z_t x_t - x_signal # 4. 噪声分析 # (a) 均值和方差 z_mean np.mean(z_t) z_variance np.var(z_t) # (b) 自相关函数 def autocorrelation(x): # 使用 full 模式得到完整的自相关函数 corr correlate(x, x, modefull, methodauto) # 截取正时延部分自相关函数在负时延部分是对称的 return corr[len(x)-1:] z_autocorr autocorrelation(z_t) # 对自相关函数进行归一化使得零延迟处的自相关值为 1 z_autocorr_normalized z_autocorr / z_autocorr[0] # (c) 功率谱密度 N len(z_t) dt time[1] - time[0] # 采样间隔 frequencies fftfreq(N, dt) z_fft fft(z_t) z_psd np.abs(z_fft)**2 / N # 标准化功率谱密度 # 5. 可视化 # 设置图形显示 fig, axs plt.subplots(3, 2, figsize(14, 12)) # (1) 绘制接收到的信号 x(t) 和已知信号 x_signal axs[0, 0].plot(time, x_t, labelReceived Signal x(t), colorblue) axs[0, 0].plot(time, x_signal, labelKnown Signal x_signal, linestyle--, colororange) axs[0, 0].set_title(Received Signal and Known Signal) axs[0, 0].set_xlabel(Time [s]) axs[0, 0].set_ylabel(Signal Amplitude) axs[0, 0].legend() # (2) 绘制噪声部分 z(t) axs[0, 1].plot(time, z_t, labelNoise z(t), colorred) axs[0, 1].set_title(Noise z(t)) axs[0, 1].set_xlabel(Time [s]) axs[0, 1].set_ylabel(Noise Amplitude) # (3) 绘制噪声自相关函数归一化后的自相关函数 axs[1, 0].plot(np.arange(0, len(z_autocorr_normalized)), z_autocorr_normalized, colorgreen) axs[1, 0].set_title(Autocorrelation of Noise (Normalized)) axs[1, 0].set_xlabel(Lag) axs[1, 0].set_ylabel(Autocorrelation) # (4) 绘制噪声的功率谱密度 axs[1, 1].plot(frequencies[:N // 2], z_psd[:N // 2], colorpurple) axs[1, 1].set_title(Power Spectral Density of Noise) axs[1, 1].set_xlabel(Frequency [Hz]) axs[1, 1].set_ylabel(Power) # (5) 绘制接收到的信号 x(t) 和噪声部分 z(t) axs[2, 0].plot(time, x_t, labelReceived Signal x(t), colorblue) axs[2, 0].plot(time, z_t, labelNoise z(t), colorred) axs[2, 0].set_title(Received Signal and Noise z(t)) axs[2, 0].set_xlabel(Time [s]) axs[2, 0].set_ylabel(Amplitude) axs[2, 0].legend() # (6) 显示统计量均值和方差 axs[2, 1].axis(off) text fNoise Mean: {z_mean:.3e}\nNoise Variance: {z_variance:.3e} axs[2, 1].text(0.1, 0.5, text, fontsize12) # 布局调整 plt.tight_layout() plt.show()问题二import numpy as np import pandas as pd import matplotlib.pyplot as plt from scipy.signal import butter, filtfilt # 读取Excel文件数据 file_path rD:\A题\pythonProject1\第二题.xlsx # Excel文件路径 data pd.read_excel(file_path) # 提取时间和信号列 time data[Received Signal Time].values signal data[Received Signal Value].values # 确定采样时间间隔 Ts Ts time[1] - time[0] # 信号去噪使用低通滤波器 def lowpass_filter(data, cutoff_frequency, sample_rate, order4): nyquist 0.5 * sample_rate normal_cutoff cutoff_frequency / nyquist b, a butter(order, normal_cutoff, btypelow, analogFalse) filtered_data filtfilt(b, a, data) return filtered_data # 假设信号的频率大约在1e7 Hz左右可以选择一个适当的截止频率 cutoff_freq 1e7 # 截止频率单位Hz filtered_signal lowpass_filter(signal, cutoff_freq, 1 / Ts) # 计算信号的傅里叶变换 N len(filtered_signal) # 数据点的个数 frequencies np.fft.fftfreq(N, Ts) # 获取频率对应的数组 fft_values np.fft.fft(filtered_signal) # 计算傅里叶变换 # 取傅里叶变换的模得到信号的频谱 fft_magnitude np.abs(fft_values) # 去掉负频率部分 positive_frequencies frequencies[:N // 2] positive_fft_magnitude fft_magnitude[:N // 2] # 对频谱的幅度使用对数缩放以便更清晰地显示 log_fft_magnitude np.log10(positive_fft_magnitude 1e-10) # 1e-10 防止对数为负无穷 # 找到频谱中的最大峰值对应的频率 peak_frequency_index np.argmax(positive_fft_magnitude) estimated_frequency positive_frequencies[peak_frequency_index] # 输出估计的频率 print(f估计的信号频率为{estimated_frequency:.2e} Hz) # 可视化时域信号和频谱 plt.figure(figsize(12, 8)) # 时域信号图 plt.subplot(3, 1, 1) plt.plot(time, signal, labelOriginal Signal) plt.plot(time, filtered_signal, labelFiltered Signal, linestyle--, colororange) plt.xlabel(Time (s)) plt.ylabel(Amplitude) plt.title(Received Signal in Time Domain) plt.legend() plt.grid(True) # 频谱图线性幅度 plt.subplot(3, 1, 2) plt.plot(positive_frequencies, positive_fft_magnitude, labelFFT Magnitude, colororange) plt.xlabel(Frequency (Hz)) plt.ylabel(Magnitude) plt.title(Frequency Spectrum (Linear Scale)) plt.axvline(xestimated_frequency, colorr, linestyle--, labelfEstimated Frequency: {estimated_frequency:.2e} Hz) plt.legend() plt.grid(True) # 频谱图对数幅度 plt.subplot(3, 1, 3) plt.plot(positive_frequencies, log_fft_magnitude, labelLog of FFT Magnitude, colorpurple) plt.xlabel(Frequency (Hz)) plt.ylabel(Log Magnitude) plt.title(Frequency Spectrum (Log Scale)) plt.axvline(xestimated_frequency, colorr, linestyle--, labelfEstimated Frequency: {estimated_frequency:.2e} Hz) plt.legend() plt.grid(True) plt.tight_layout() plt.show()问题三import numpy as np import matplotlib.pyplot as plt import pandas as pd # 1. 从Excel文件中读取数据 file_path rD:\A题\pythonProject1\第三题.xlsx data pd.read_excel(file_path) # 2. 提取时间和信号值数据 time data[Received Signal Time].values signal data[Received Signal Value].values # 3. 获取采样间隔 T_s time[1] - time[0] # 计算时间间隔 # 4. 进行零填充使得FFT的分辨率更高 N len(signal) # 信号长度 N_padded 2**np.ceil(np.log2(N)).astype(int) # 选择大于等于N的最小2的幂次增加FFT的分辨率 signal_padded np.pad(signal, (0, N_padded - N), constant, constant_values(0, 0)) # 5. 对信号进行快速傅里叶变换FFT fft_signal np.fft.fft(signal_padded) # 计算FFT frequencies np.fft.fftfreq(N_padded, T_s) # 生成频率轴 # 6. 计算频谱的幅度只考虑正频率部分 magnitude np.abs(fft_signal[:N_padded // 2]) frequencies frequencies[:N_padded // 2] # 只取正频率部分 # 7. 对频谱进行平滑处理减小噪声对频谱的影响 def smooth_spectrum(magnitude, window_size5): 使用滑动平均法对频谱进行平滑处理 return np.convolve(magnitude, np.ones(window_size)/window_size, modesame) # 平滑后的频谱 smoothed_magnitude smooth_spectrum(magnitude) # 8. 使用对数尺度来展示频域数据对于幅度的对数 log_magnitude np.log(smoothed_magnitude 1e-10) # 防止出现log(0)的情况 # 9. 找到主峰对应的频率 peak_frequency frequencies[np.argmax(smoothed_magnitude)] # 找到幅度最大值对应的频率 # 输出估计的频率 print(f估计的信号频率{peak_frequency:.2e} Hz) # 10. 可视化时域信号和频域信号 fig, axs plt.subplots(2, 1, figsize(10, 8)) # 时域信号 axs[0].plot(time, signal, labelReceived Signal, colorb) axs[0].set_title(Received Signal in Time Domain) axs[0].set_xlabel(Time (s)) axs[0].set_ylabel(Signal Value) axs[0].grid(True) # 频域信号使用对数尺度 axs[1].plot(frequencies, log_magnitude, labelLog Magnitude of Smoothed FFT, colorr) axs[1].set_title(Signal Frequency Spectrum (Log Scale)) axs[1].set_xlabel(Frequency (Hz)) axs[1].set_ylabel(Log Magnitude) axs[1].grid(True) # 显示估计的频率 axs[1].axvline(xpeak_frequency, colorg, linestyle--, labelfEstimated Frequency: {peak_frequency:.2e} Hz) axs[1].legend() # 显示图形 plt.tight_layout() plt.show()问题四import numpy as np import pandas as pd import matplotlib.pyplot as plt # 1. 从 Excel 文件中读取数据 file_path rD:\A题\pythonProject1\第四题.xlsx # 假设 Excel 文件的第一张工作表Sheet1包含我们需要的数据 # 数据列假设为 Received Signal Time 和 Received Signal Value data pd.read_excel(file_path, sheet_name0) # 提取时间和信号值 time_data data[Received Signal Time].to_numpy() signal_data data[Received Signal Value].to_numpy() # 打印前几行数据以确认读取正确 print(读取的数据) print(data.head()) # 2. 间歇接收分析 # 计算时间差 time_diff np.diff(time_data) # 时间间隔 avg_time_diff np.mean(time_diff) # 平均时间间隔 # 打印时间间隔信息 print(f时间间隔{time_diff}) print(f平均时间间隔{avg_time_diff}) # 3. 使用FFT估计信号频率 # 我们将信号数据进行傅里叶变换以估计信号的频率成分 sampling_rate 1 / avg_time_diff # 采样频率基于时间间隔的倒数 n len(signal_data) # 信号的长度 frequency np.fft.fftfreq(n, davg_time_diff) # 计算频率轴 fft_signal np.fft.fft(signal_data) # 傅里叶变换 # 计算功率谱频域的幅度 power_spectrum np.abs(fft_signal) ** 2 # 4. 可视化结果 # 时域信号图 plt.figure(figsize(12, 6)) plt.subplot(2, 1, 1) plt.plot(time_data, signal_data, labelReceived Signal, colorb) plt.title(Time Domain: Received Signal) plt.xlabel(Time (s)) plt.ylabel(Signal Amplitude) plt.grid(True) # 频域信号图 plt.subplot(2, 1, 2) plt.plot(frequency[:n//2], power_spectrum[:n//2], labelPower Spectrum, colorr) plt.title(Frequency Domain: Power Spectrum) plt.xlabel(Frequency (Hz)) plt.ylabel(Power) plt.grid(True) plt.tight_layout() plt.show() # 5. 频率估计 # 频率估计可以通过找到功率谱的峰值频率来实现 peak_frequency_idx np.argmax(power_spectrum[:n//2]) # 获取最大功率谱的索引 estimated_frequency abs(frequency[peak_frequency_idx]) # 对应的频率值 print(f估计的频率{estimated_frequency} Hz)