基于低通滤波器的原子钟时差测量数据降噪处理

1. Savitzky-Golay 滤波器频率响应图，采用窗口长度为 51、多项式阶数为 3 的参数设置

# Generate filter coefficients for an impulse response
impulse = np.zeros(window_length)
impulse[window_length // 2] = 1  # Set a single impulse at the center
h = savgol_filter(impulse, window_length, polyorder)  # Generate filter coefficients

# Compute frequency response
w, H = freqz(h, worN=8000)

# Plot frequency response
plt.figure(figsize=(10, 6))
plt.plot(w / np.pi, 20 * np.log10(np.abs(H)), label='Magnitude Response')
plt.title("Savitzky-Golay Filter Frequency Response (Window=51, Polyorder=3)")
plt.xlabel("Normalized Frequency (×π rad/sample)")
plt.ylabel("Magnitude (dB)")
plt.ylim(-100, 5)
plt.grid(True)
plt.legend()
plt.show()

2. 归一化截止频率 fc 随 M变化的趋势图

import numpy as np
import matplotlib.pyplot as plt

# Parameters
M1 = np.linspace(25, 100, 100)  # M >= 25
M2 = np.linspace(10, 24, 100)   # 10 <= M < 25
N = 20                          # Example N value, N < M

# Compute fc for both cases
fc1 = (N + 1) / (3.2 * M1 - 4.6)  # For M >= 25
fc2 = (N + 1) / (3.2 * M2 - 2)    # For 10 <= M < 25

# Plot
plt.figure(figsize=(10, 6))
plt.plot(M1, fc1, label=r'$M \geq 25$', color='blue', linewidth=2)
plt.plot(M2, fc2, label=r'$10 \leq M < 25$', color='orange', linewidth=2)
plt.title("Variation of Normalized Cutoff Frequency $f_c$ with $M$", fontsize=14)
plt.xlabel("$M$", fontsize=12)
plt.ylabel("$f_c$", fontsize=12)
plt.axvline(25, color='gray', linestyle='--', label='$M=25$ Boundary')
plt.grid(True)
plt.legend(fontsize=12)
plt.show()

3. 模拟实验数据的时间差变化图

# Parameters
sampling_interval = 1  # Sampling interval in seconds
time_span = 4 * 24 * 3600  # Time span of 4 days in seconds
time = np.arange(0, time_span, sampling_interval)  # Time array

# Simulated phase difference data (random walk noise for phase differences)
np.random.seed(42)
phase_diff = np.cumsum(np.random.normal(0, 0.1, size=len(time)))  # Simulated random walk for phase

# Simulated time difference data with uncertainty
rise_time_slow = 500e-12  # 500 ps uncertainty for slower rise time pulses
rise_time_fast = 100e-12  # 100 ps uncertainty for faster rise time pulses
time_diff = phase_diff + np.random.normal(0, rise_time_fast, size=len(time))  # Adding fast pulse noise

# Introduce slower rise time regions (higher uncertainty)
slow_rise_region = (time > 1 * 24 * 3600) & (time < 2 * 24 * 3600)  # Second day
time_diff[slow_rise_region] += np.random.normal(0, rise_time_slow, size=slow_rise_region.sum())

# Simulate additional uncertainty due to phase comparator (5.1 × 10^-13 added uncertainty)
additional_uncertainty = 5.1e-13  # Uncertainty due to phase comparator
time_diff += np.random.normal(0, additional_uncertainty, size=len(time))

# Plot the data
plt.figure(figsize=(12, 6))
plt.plot(time / 3600, time_diff, label="Time Difference Data", color="blue", linewidth=1)
plt.title("Simulated Time Difference Data with Phase Comparator Effects", fontsize=14)
plt.xlabel("Time (hours)", fontsize=12)
plt.ylabel("Time Difference (s)", fontsize=12)
plt.axvspan(24, 48, color='orange', alpha=0.3, label="Higher Uncertainty Region (Slow Pulse Rise)")
plt.grid(True)
plt.legend(fontsize=12)
plt.tight_layout()
plt.show()

4. 相差数据与时差数据模拟图

# Parameters for simulation
sampling_interval = 1  # Sampling interval in seconds
time_span = 4 * 24 * 3600  # 4 days in seconds
time = np.arange(0, time_span, sampling_interval)  # Time array in seconds

# Convert time to reduced Julian date (RJD)
rjd_start = 59292  # Start of the reduced Julian date
rjd_time = rjd_start + time / (24 * 3600)  # Convert seconds to RJD

# Simulate phase difference data (linear trend with random noise)
np.random.seed(42)
phase_diff = -1e7 * (rjd_time - rjd_start) + np.random.normal(0, 0.1, size=len(time))

# Simulate time difference data (based on phase difference with added noise)
time_diff = phase_diff + np.random.normal(0, 0.05, size=len(time))  # Smaller noise for time diff

# Plot full data: phase difference and time difference
plt.figure(figsize=(12, 6))
plt.plot(rjd_time, phase_diff, label="Phase Difference Data", color="red", linewidth=1)
plt.plot(rjd_time, time_diff, label="Time Difference Data", color="blue", linewidth=1)
plt.title("Phase Difference and Time Difference Measurement Results", fontsize=14)
plt.xlabel("Reduced Julian Date (RJD)", fontsize=12)
plt.ylabel("Measurement Result (×10⁷) s", fontsize=12)
plt.legend(fontsize=12)
plt.grid(True)
plt.tight_layout()
plt.show()

# Plot zoomed-in data
zoom_region = (rjd_time >= 59293) & (rjd_time <= 59293.5)  # Define zoom region
plt.figure(figsize=(12, 6))
plt.plot(rjd_time[zoom_region], phase_diff[zoom_region], label="Phase Difference Data", color="red", linewidth=1)
plt.plot(rjd_time[zoom_region], time_diff[zoom_region], label="Time Difference Data", color="blue", linewidth=1)
plt.title("Zoomed-in View of Phase Difference and Time Difference Data", fontsize=14)
plt.xlabel("Reduced Julian Date (RJD)", fontsize=12)
plt.ylabel("Measurement Result (×10⁷) s", fontsize=12)
plt.legend(fontsize=12)
plt.grid(True)
plt.tight_layout()
plt.show()

5. 不同多项式阶数与交叉验证误差图

# 模拟不同多项式阶数与交叉验证误差数据
window_widths = np.arange(60, 360, 10)  # 窗口宽度范围
poly_degrees = [1, 2, 3, 4]  # 多项式阶数

# 生成模拟交叉验证误差数据
errors = {}
np.random.seed(42)  # 固定随机种子以确保结果一致性
for degree in poly_degrees:
    errors[degree] = 0.6 + np.sin(window_widths / (50 + degree * 20)) * 0.005 + degree * 0.0005

# 绘制不同多项式阶数下的交叉验证误差图
plt.figure(figsize=(12, 8))
for degree in poly_degrees:
    plt.plot(window_widths, errors[degree], label=f"Polynomial Degree {degree}", linewidth=1.5)

plt.title("Cross-validation Error vs Window Width for Different Polynomial Degrees", fontsize=14)
plt.xlabel("Window Width", fontsize=12)
plt.ylabel("Cross-validation Error", fontsize=12)
plt.legend(fontsize=10, loc="upper right")
plt.grid(True, linestyle="--", linewidth=0.5)
plt.tight_layout()
plt.show()