0

出于学习目的，我一直在尝试在 Python 中实现 Prof.Selesnick 的 DWT 实现的 MATLAB 代码。

function [lo, hi] = afb(x, af)

% [lo, hi] = afb(x, af)
%
% Analysis filter bank
% x -- N-point vector (N even); the resolution should be 2x filter length
%
% af  -- analysis filters
% af(:, 1): lowpass filter (even length)
% af(:, 2): highpass filter (even length)
%
% lo: Low frequency
% hi: High frequency
%

N = length(x);
L = length(af)/2;
x = cshift(x,-L);

% lowpass filter
lo = upfirdn(x, af(:,1), 1, 2);
lo(1:L) = lo(N/2+[1:L]) + lo(1:L);
lo = lo(1:N/2);

% highpass filter
hi = upfirdn(x, af(:,2), 1, 2);
hi(1:L) = hi(N/2+[1:L]) + hi(1:L);
hi = hi(1:N/2);

我特别困在lo(1:L) = lo(N/2+[1:L]) + lo(1:L);

我尝试lo[np.arange(0,L)]=lo[N // 2 + np.concatenate([np.arange(0,L)])) + lo[np.arange(0,L)] 过，但似乎没有用。将不胜感激任何帮助。

我有一个大小为 10,000 的输入信号 x，当我执行代码时，它会停在该特定行并显示index 5001 is out of bounds for axis 0 with size 5001. 我好像越界了

4

3 回答 3

0

在 Octave 会话中：

>> lo = 1:10;
>> L=3; N=4;
>> N/2+[1:L]
ans =

   3   4   5

>> lo(N/2+[1:L])+lo(1:L)
ans =

   4   6   8

在numpy：

In [100]: lo = np.arange(1,11)
In [101]: L=3; N=4
In [102]: N/2+np.arange(0,L)
Out[102]: array([2., 3., 4.])
In [105]: lo[int(N/2)+np.arange(0,L)]+lo[:L]
Out[105]: array([4, 6, 8])

等效地

In [106]: n=int(N/2); lo[n:n+L]+lo[:L]
Out[106]: array([4, 6, 8])
于 2021-08-30T04:55:22.920 回答
0

numpy 库已经有一堆可以帮助我们处理切片的项目。

例如，当您需要对 numpy 数组进行a = np.array([1,2,3,4,5,6,7,8,9,10])切片时，您可以使用切片a[:n]并将n其作为序列中的最后一个元素。如果n = 2那么a[:2] = [1,2]. 试试这个，你会写出干净的代码。

对于您的“越界”错误，我建议您检查索引并记住 Python 计数从 0 开始。

例如（使用我们上面定义的“a”数组）：

a[0] = 1
a[-1] = 10 
a[9] = 10

请注意，最后一个元素位于索引9处。我相信这是您的错误，您必须使用foo[5000]而不是foo[5001].

我希望我对你有一点帮助！

于 2021-08-30T01:30:50.740 回答
0

我不熟悉该upfirdn功能，但看起来您的功能应用了前向 DWT 的一级？

这个答案包含 Python 中（多级）DWT 的前向实现。这是 Mallat 的古老 FWT 配方，可能与您的略有不同，但基于相同的原理：

  • 前向 FWT：级别i-> 过滤器 -> 下采样 -> 级别i+1
  • 逆 FWT：级别i+1-> 上采样 -> 过滤器 -> 级别i
于 2021-08-30T13:11:46.400 回答