数学算法所有主题 HTML Java 中级开发 Windows Linux

离散小波变换，一种 Java 实现

Mark H Bishop

4.95/5 (16投票s)

2014年10月4日

CPOL

10分钟阅读

62831

5712

本文介绍了一个执行离散小波变换的 Java 示例应用程序。

DwtSrc.zip sha256 和: 665bfbd806214149bd569571255e1cc15ffd660124cede54caf22bf7964db86c

DwtRunnable.zip sha256 和: 95030fe0f023bcefc9874c5580b42dddb6acfe0f14bc9caf3d6c635f5b844bd4

引言

小波变换是理解世界（又称信号处理）的重要工具。术语信号指任何随某个范围变量函数变化的信息流。通常，在信号处理文献中，范围变量是时间，但它本质上可以是任何东西：位置、施加的磁场、波长，甚至是无量纲序列。在应用中，该函数可以是记录的声音、图像、传感器响应、心电图或核磁共振、光谱数据、气候记录、混沌等。这个列表是无穷无尽的，开发者有无数的创意机会来探索使用小波进行数据分析的新方法。一个非常活跃的学科是基于机器的模式识别。

在本文中，我提供了一个使用离散小波变换来探索一维信号的应用程序。该方法可以进一步开发以变换更高维度的信号，如图像。还有其他方法可以将小波功能添加到 Java 中，例如使用 Matlab/Scilab-Wavelab 包装器或开源库。为了说明，我选择直接在 Java 方法中展示计算的基础。

本项目中的计算类是我几年前编写的 C^# 代码的翻译。反向过程应该相当容易。

背景

您可能是在 Google 上搜索类似“计算小波变换”后才看到本文的。一些搜索结果提供了重要的定义和不同复杂程度的数学背景。另一些则通过简单地说明“打开 Matlab/Scilab/等并输入...”来回答如何计算变换的问题。

这些答案是有效的，但开发者希望看到或找出他们选择的语言中算法的实现。以下是我的实现。有更优雅和高效的方法，但我选择字面地在单个 CPU 线程上运行的方法中表达算法。

开发此演示的一个挑战是多分辨率数据的有效可视化。打印演示需要使用堆叠图，这在计算机显示器上进行交互式评估是很困难的。在这种情况下，我选择在同一张图上同时绘制刻度，允许用户从菜单项中显示/隐藏单个刻度。此外，我将数据结构化，使得所有刻度数据的叠加将完美地重建信号。

尝试可运行的 jar

要运行演示，请下载上面的 DwtRunnable.zip，解压它，找到 DwtRunnable.jar 文件并运行它。在 Windows 上，我只需双击 DwtRunnable.jar 文件即可运行它。在 Linux 上，我在终端中输入 java -jar <thePath>/DwtRunnable.jar。显然，您需要安装 Java，最好是最新版本并更新所有补丁。您可以从主 GUI 运行示例。单击 Plot 查看图表和信号的数据表。从数据表窗口，您可以执行变换或多分辨率分析。您可以加载菜单示例，也可以在表达式框中输入一个函数。此外，如果单击 Empty Table，将出现一个空白表，您可以在其中打开一个 csv 文件进行绘图和分析。

以下是使用高斯修改混合正弦示例的一些屏幕截图。

这是初始 GUI、其图表和其数据表

Entry GUI Function plot Data table

这是变换输入对话框、离散小波变换及其逆变换（重建）

Parameter dialog Discrete transform Inverse transform

最后，这是多分辨率分析及其表格

Co-plot of MRA MRA table

多分辨率表中的列将总和为原始信号。MRA 图中的各个刻度可以从菜单中显示/隐藏。

使用代码

DwtSrc.zip 下载是一个压缩的 Eclipse 项目文件夹。如果您使用 Eclipse，请将其解压到您选择的工作区中，并从 IDE 中创建一个名为 DwtDemo 的项目。这应该会填充项目并启用测试运行。对于其他 IDE，或者如果您选择其他方式创建 Eclipse 项目，您将需要将 jel.jar 库链接到您的构建路径中。它位于 jel-2.0.1 文件夹中。有一篇关于使用 JEL 表达式库的先前文章，位于 https://codeproject.org.cn/Articles/737497/Java-Function-Compiler。它是演示的一部分，但对于使用离散小波变换特定类来说不是必需的。

包括以下 16 个类

枚举

Colors.java（一些能很好地绘制多分辨率图的颜色）
Wavelet.java（Haar、Daubechies、Symmlet、Coiflet 等）

基于 JFrame 的交互式 GUI

MainForm.java（演示入口点）
MatrixForm.java（数据表和相关菜单）
ParameterInput.java（用户选择变换参数）
PlotFrame.java（带菜单的绘图）

计算类

DWT.java（文章主要主题）
OrthogonalFilters.java（也主要）
MatrixOps.java（由 DWT.java 使用）

其他实用程序

VisualizeMultiResolution.java（应用程序特定帮助程序）
Examples.java（一些特殊函数，一些非解析函数）
FileOps.java
ImagePanel.java
Signal.java（用于管理 xy 数据的应用程序特定帮助程序）
StringUtils.java
VariableProvider.java（JEL 需要）

您会在每个类开始的地方找到一个注释，说明该类的职责。我尝试在需要的地方对方法进行适当的注释。错误处理是最小的。关于 DWT.java 和 OrthogonalFilters.java，有几件重要的事情需要讨论。

DWT.java 包含方法 transform()

public static double[] transform(double[] signal, Wavelet wavelet,
            int order, int L, Direction direction) throws Exception {
        if (direction == Direction.forward) {
            return forwardDwt(signal, wavelet, order, L);
        } else {
            return inverseDwt(signal, wavelet, order, L);
        }
    }

信号是您想要变换的序列，其长度必须是 2 的偶数幂。（如果需要，还提供了另一个用于填充的方法。）
小波可以是 Haar、Daubechies 等。
阶数是小波特定的选择，必须与 OrthogonalFilters.java 中为每种小波类型定义的选择之一相符。这些滤波器名称的命名法和形式在文献中差异很大。
L 是您想要包含在变换中的最粗糙的尺度。
方向对于正向变换是forward，对于逆变换是reverse（在类中定义）。

根据方向，它将调用 forwardDwt() 或 inverseDwt()。

对于正向变换，我更喜欢迭代地为每个尺度形成一个正交镜像滤波器（QMF），使用 OrthogonalFilters.java 中提供的适当高通和低通滤波器作为矩阵。镜像矩阵在每个尺度上由 makeQMFMatrix() 创建，级联的矩阵-向量乘法实现变换。返回结果 dWT 最初设置为输入信号。每次尺度减半时，包含 dWT 前半部分的向量（subResult）乘以适当大小的镜像矩阵。随后，乘积覆盖 dWT 的前 n/(2ⁱ) 个值，其中 i 是 for 循环中的当前迭代。

		public static double[] forwardDwt(double[] signal, Wavelet wavelet,
            int order, int L) throws Exception {
        int n = signal.length;
        if (!isValidChoices(wavelet, order, L, n)) {
            throw new Exception(
                    "Invalid wavelet /order/scale/signal-length combination.");
        }
        double[] dWT = MatrixOps.deepCopy(signal);
        int log2n = (int) (Math.log(n) / Math.log(2));
        int iterations = log2n - L;
        int subLength = n;
        double[] H = OrthogonalFilters.getLowPass(wavelet, order);
        double[] G = OrthogonalFilters.getHighPass(H);
        for (int i = 0; i < iterations; i++) {
            subLength = n / (int) (Math.pow(2, i));
            double[][] QMF = makeQMFMatrix(subLength, H, G);
            double[] subResult = new double[subLength];
            subResult = subCopy(dWT, subResult, subLength);
            double[] temp = MatrixOps.multiply(QMF, subResult);
            dWT = subCopy(temp, dWT, subLength);
        }
        return dWT;
    }

inverseDwt() 方法可以应用于任何信号，但通常它用于重建目的，应用于先前获得的正向变换，通常在对变换进行一些修改（如去噪）之后。如果进行重建，请使用与原始变换中使用的相同的小波、阶数和最粗尺度（L）。该方法以与 forwardDwt() 类似的方式反转正向变换的操作，只是在 for (i = L + 1; i <= log₂n, i++) 循环结构中使用了 QMF 矩阵的转置作为乘数。

	public static double[] inverseDwt(double[] signal, Wavelet wavelet,
            int order, int L) throws Exception {
        int n = signal.length;
        if (!isValidChoices(wavelet, order, L, n)) {
            throw new Exception(
                    "Invalid wavelet /order/scale/signal-length combination.");
        }
        int log2n = (int) (Math.log(n) / Math.log(2));
        int subLength;
        double[] preserveCopy = new double[signal.length];
        preserveCopy = subCopy(signal, preserveCopy, signal.length);
        double[] H = OrthogonalFilters.getLowPass(wavelet, order);
        double[] G = OrthogonalFilters.getHighPass(H);
        for (int i = L + 1; i <= log2n; i++) {
            subLength = (int) (Math.pow(2, i));
            double[][] QMF = makeQMFMatrix(subLength, H, G);
            QMF = MatrixOps.transpose(QMF);
            double[] subResult = new double[subLength];
            subCopy(signal, subResult, subLength);
            subResult = MatrixOps.multiply(QMF, subResult);
            signal = subCopy(subResult, signal, subLength);
        }
        double[] iDWT = new double[n];
        iDWT = subCopy(signal, iDWT, n);
        signal = preserveCopy;
        return iDWT;
    }

有许多地方可以示意性地回顾上面实现的底层算法，包括：http://users.rowan.edu/~polikar/WAVELETS/WTpart4.html 图 4.1。这里使用的 QMF 矩阵实现更难找到，所以请自行确信它是功能正确的（下面提供了一些测试）。演示的实现效率有点低，因为它在实现时显式地将稀疏矩阵当作密集矩阵来处理。概念模型，使用矩阵-向量乘法进行滤波，非常直观，可以使用稀疏乘法策略重新实现，可能在并行 GPU 线程或 GPU 上。

多分辨率分析

执行多分辨率分析的方法 mRA() 接受与 transform() 方法相同的参数集，并返回一个 ArrayList<Object>，其规范如下

mRA.get(0) 是一个 ArrayList<double[]>，其中包含从最精细尺度到最粗糙尺度的多分辨率尺度数据，以及包含对应于大于 L 的尺度的曲线的最终近似尺度。
mRA.get(1) 是一个 int[]，包含值 j，使得 MRA 中从最精细到最粗糙的每个尺度都可以稍后使用公式：scale = 2^-j 来确定
mRA.get(0) 中列出的最后一个 double[] 和 mRA.get(1) 提供的最后一个 j 对应于近似曲线，它承载的频率大于 2^-L

下面是 mRA() 方法。最初，信号被正向变换为 dwt。第一个循环通过对 dwt 中连续的二进尺度进行零插值，然后逆变换中间结果，创建名为 mRA 的 ArrayList<double[]> 中的前 n-1 个项目。第二个循环添加了频率大于 2^-L 的最终近似曲线。第三个循环将 j 值添加到名为 scalesUsed 的 int[] 中，将最终近似曲线指定为 j = 0。最后，将 mRA 和 scalesUsed 作为对象添加到名为 result 的返回值中。

		public static ArrayList mRA(double[] signal, Wavelet wavelet,
			int order, int L) throws Exception {
		ArrayList result = new ArrayList();
		int n = signal.length;
		if (!isValidChoices(wavelet, order, L, n)) {
			throw new Exception(
					"Invalid wavelet /order/scale/signal-length combination.");
		}
		int J = (int) (Math.log(n) / Math.log(2));
		double[] dwt = forwardDwt(signal, wavelet, order, L);
		ArrayList mRA = new ArrayList();
		for (int j = (J - 1); j >= L; j--) {
			double[] w = new double[n];
			int[] dyad = dyad(j);
			for (int k = dyad[0]; k <= dyad[dyad.length - 1]; k++) {
				w[k - 1] = dwt[k - 1];
			}
			mRA.add(inverseDwt(w, wavelet, order, L));
		}
		// All frequencies lower than those revealed at L
		double[] w = new double[n];
		int limit = (int) Math.pow(2, L);
		for (int i = 0; i < limit; i++) {
			w[i] = dwt[i];
		}
		mRA.add(inverseDwt(w, wavelet, order, L));

		int[] scalesUsed = new int[mRA.size()];
		int scaleCounter = 0;
		for (int j = (J - 1); j >= L; j--) {
			int[] dyad = DWT.dyad(j);
			scalesUsed[scaleCounter] = (int) (Math.log(dyad.length) / Math
					.log(2));
			scaleCounter++;
		}
		// Next line: 0 is a dummy value for the approximation carrying all
		// lower frequencies corresponding to scales larger than 2^-L
		scalesUsed[scaleCounter] = 0;
		result.add(mRA);
		result.add(scalesUsed);
		return result;
	}

代码测试

我以几种方式测试了我的结果。

首先，我使用各种小波和小波参数变换了几个函数，检查结果以确保其符合预期，然后对结果进行逆变换并验证完美重建。在 DemoSrc 下载目录中有一个名为 TestSignal.xls 的电子表格，它呈现了一个三频混合正弦函数，在 x = 0.0525 和 x = 0.0625 之间有瞬态噪声，并在 0.154 和 0.155 处有两个尖峰。此信号用于创建名为 TestSignalXY.csv 的 csv 文件。

Test signal

DWT.java 中有一个名为 relativeErrors() 的测试方法，它能够自动对 OrthogonalFilters.java 中提供的所有正交滤波器组合进行重建测试。使用 TestSignalXY.csv 作为数据文件运行此方法，会显示原始信号与每个小波/阶数组合的重建信号之间差异的相对残差。除了 Battle 小波约为 ~10^-4 之外，所有小波的值都在 ~10^-12 到 ~10^-16 之间。Battle 小波的系数仅提供到 9 位有效数字，因此重建效果预计会降低。使用该应用程序，MRA 矩阵的列相加可以完美地重现信号。下面以堆叠格式显示了选定尺度的多分辨率图。最精细的尺度显示尖峰，中等尺度显示瞬态，其余图显示三个正弦波。

Stacked MRA