我正在尝试编写代码来实现离散小波变换（haar 小波 dwt），而不使用 python 中的包。
到目前为止，我找到了一个链接，他们实现了类似的东西，链接这个小波变换实现正确吗？。运行时没有报错，但最终结果不正确。我运行的代码是：
def discreteHaarWaveletTransform(x):
N = len(x)
output = [0.0]*N
length = N >> 1
while True:
for i in range(0,length):
summ = x[i * 2] + x[i * 2 + 1]
difference = x[i * 2] - x[i * 2 + 1]
output[i] = summ
output[length + i] = difference
if length == 1:
return output
#Swap arrays to do next iteration
x = output[:length << 1]
length >>= 1
输入：
list=[56, 40, 8, 24, 48, 48, 40, 16]
电流输出：
[280, -24, 64, 40, 16, -16, 0, 24]
预期输出：
[35, -3, 16, 10, 8, -8, 0, 12]
有什么明显的东西我看不到吗？
类似这样的东西应该可以解决问题——它几乎是这个 C 答案的直译。这可能意味着这不是非常 Python-y 的代码，但如果你没有使用 numpy
来实现它，那无论如何也不是非常 Python-y。您没有获得预期输出的主要原因是您忘记在过滤后缩放输出。这使得每个下一级的系数大约高两倍。
请注意，1/2 的缩放可以提供预期的输出，但更常用的是 1/2√2 的缩放，以保留小波变换下信号的 L2 范数。
def haarFWT ( signal, level ):
s = .5; # scaling -- try 1 or ( .5 ** .5 )
h = [ 1, 1 ]; # lowpass filter
g = [ 1, -1 ]; # highpass filter
f = len ( h ); # length of the filter
t = signal; # 'workspace' array
l = len ( t ); # length of the current signal
y = [0] * l; # initialise output
t = t + [ 0, 0 ]; # padding for the workspace
for i in range ( level ):
y [ 0:l ] = [0] * l; # initialise the next level
l2 = l // 2; # half approximation, half detail
for j in range ( l2 ):
for k in range ( f ):
y [j] += t [ 2*j + k ] * h [ k ] * s;
y [j+l2] += t [ 2*j + k ] * g [ k ] * s;
l = l2; # continue with the approximation
t [ 0:l ] = y [ 0:l ] ;
return y
def main():
s0 = [ 56, 40, 8, 24, 48, 48, 40, 16 ];
print( "level 0" );
print( s0 );
print( "level 1" );
print( haarFWT (s0, 1 ) );
print( "level 2" );
print( haarFWT (s0, 2 ) );
print( "level 3" );
print( haarFWT (s0, 3 ) );
if __name__ == "__main__":
main()
# run with: >>> execfile ( "haarwavelet.py" )
pywt.wavedec(signal, "haar", mode="zero")
的pywavelet实现一致。
def haar_wavedec(signal, level=None, scale=np.sqrt(0.5)):
"""
A numpy-based implementation of a multi-level 1D wavelet decomposition
using the Haar wavelet. The input signal must be at least 2 samples long,
and is always zero-padded to a length of (2^(level+1)) samples.
Parameters
----------
signal : np.ndarray
A numpy array containing the input signal.
level : int or None, optional
A number of levels of the decomposition or
None for automatic estimation.
The default is None.
scale : float, optional
A scaling factor for the decomposition.
The default value preserves L2 norm of the transformed signal.
The default is np.sqrt(0.5).
Returns
-------
res : List[np.ndarray]
A list of numpy arrays containing the result in pywavelets format:
[cA(n), cD(n), cD(n-1), ..., cD(1)], i.e., the approximate coefficients
from the n-th level and all the detail coefficients ordered
from the n-th level to the 1st.
"""
if len(signal) < 2: return None
# get the level n (by maximizing n while 2^n <= signal length) if not specified
n = int(np.log2(len(signal))) if level is None else level
h = np.array([1.0, 1.0]) # low-pass filter (sum)
g = np.array([1.0, -1.0]) # high-pass filter (difference)
# zero-pad the signal to 2^(n+1) length
t = np.pad(signal, (0,2**(n+1)-len(signal)), constant_values=0)
res = list()
for i in range(n, 0, -1):
l = len(t) // 2
# array c for intermediate results of the convolution
c = np.zeros(l * 2)
for j in range(l):
for k in range(len(h)):
c[j] += t[2*j + k] * h[k] * scale # calculate the low-pass convolution
c[j+l] += t[2*j + k] * g[k] * scale # calculate the high-pass convolution
# the next iteration will process the low-pass coefficients
t = c[:l]
# store the high-pass (detail) coefficients (stripped zeros from right)
res.append(c[l:l+np.argwhere(c[l:]).max()+1])
# store the final low-pass (approximate) coefficients
res.append(c[:l][:np.argwhere(c[:l]).max()+1])
return res[::-1]