高斯白噪声模型中Haar小波估计器实现问题求助

我尝试用Python从零实现高斯白噪声模型（参考论文公式1.1）中的Haar小波估计器，但输出结果与真实值偏差很大，怀疑存在理解误区或代码错误，同时想了解有没有现成的工具包可以直接实现该功能。

我的代码 #

import matplotlib.pyplot as plt
import numpy as np

def psi(k,l,t):
  if (l-1)/2**k<=t<l/2**k:
    return(-1)
  if  l/2**k<t<(l+1)/2**k:
    return(1)
  else:
    return(0)

def waveletcoef(k,l,X):
  S=0
  for i in range(len(X)-1):
    S=S+(X[i+1]-X[i])*psi(k,l,i/len(X))
  return(S)

def Wavelet_aprox(K,X,t):
  S=0
  for k in range(K):
    for l in range(2**k):
      S=S+waveletcoef(k,l,X)*psi(k,l,t)
  return(S)

def f(x):
  return(np.cos(x*(4*np.pi)))

m=1000
truefonc=np.zeros(m)
for i in range(m):
  truefonc[i]= f(i/m)

n=1000

x = np.zeros(m)
x[0]=f(0)
for i in range(m - 1):
    x[i + 1] = x[i] + f(i/m)*1/m + (1/(m*n))**(1/2) * np.random.randn()

K=9

approx1=np.zeros(m)
for i in range(m):
  approx1[i]= Wavelet_aprox(K,x,i/m)

plt.plot([i/m for i in range(m)],approx1, 'red')
plt.plot([i/m for i in range(m)],truefonc, 'blue')

输出结果 #

问题分析与修正 #

1. Haar小波定义错误 #

你的psi函数完全偏离了标准Haar小波的定义：

• 标准Haar小波的平移缩放形式为 $\psi_{k,l}(t) = 2^{k/2} \psi(2^k t - l)$，其中基小波$\psi(t)$在$[0,0.5)$取1，$[0.5,1)$取-1，其余为0
• 你的函数区间方向反了，且缺少归一化因子$2^{k/2}$，导致系数幅度和重构结果完全错误

修正后的小波函数：

def psi(k, l, t):
    t_scaled = 2**k * t - l
    if 0 <= t_scaled < 0.5:
        return 2**(k/2)
    elif 0.5 <= t_scaled < 1:
        return -2**(k/2)
    else:
        return 0

2. 小波系数计算的积分近似问题 #

原代码用区间左端点的小波值近似整个区间的积分，当小波在区间内跳变时会引入误差。可以改用区间中点的小波值提升精度：

def waveletcoef(k, l, X):
    dt = 1 / len(X)
    S = 0
    for i in range(len(X)-1):
        t_mid = (i + 0.5) * dt  # 取区间中点计算小波值
        S += (X[i+1] - X[i]) * psi(k, l, t_mid)
    return S

3. 重构缺少噪声阈值处理 #

根据论文模型，小波系数包含信号项和噪声项，直接重构会保留大量噪声。需要对系数做硬阈值处理过滤噪声：

def Wavelet_aprox(K, X, t):
    dt = 1 / len(X)
    sigma = (1 / (len(X)*n))**0.5
    threshold = sigma * np.sqrt(2 * np.log(len(X)))  # 通用硬阈值公式
    S = 0
    for k in range(K):
        for l in range(2**k):
            coef = waveletcoef(k, l, X)
            # 硬阈值处理：小于阈值的系数置0
            if abs(coef) > threshold:
                S += coef * psi(k, l, t)
    return S

现成工具包推荐 #

Python中最常用的小波分析库是PyWavelets（pywt），它支持Haar小波在内的多种小波变换，内置系数阈值处理和重构功能，能大幅简化实现：

• 安装：pip install pywavelets
• 核心用法示例：

import pywt

dt = 1 / m
# 对X的差分序列（近似导数）做小波变换
coeffs = pywt.wavedec(np.diff(x)/dt, 'haar', level=K)
# 计算阈值并处理系数
sigma = (1 / (m*n))**0.5
threshold = sigma * np.sqrt(2 * np.log(m))
coeffs_thresh = [pywt.threshold(c, threshold, mode='hard') for c in coeffs]
# 重构估计的f(t)
f_approx = pywt.waverec(coeffs_thresh, 'haar')

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