31 信号处理

本章介绍 Octave 中可用的信号处理和快速傅立叶变换函数。快速傅立叶变换的计算使用 FFTWFFTPACK 库,具体取决于 Octave 的构建方式。

 
y = fft (x)
y = fft (x, n)
y = fft (x, n, dim)

使用快速傅立叶变换(FFT)算法计算 x 的离散傅立叶变换。

FFT 沿着数组的第一个非单一维度进行计算。因此,如果 x 是一个矩阵,fft (x) 计算 x 每一列的 FFT。

如果使用两个参数调用,n 应为指定要使用的 x 元素数量的整数,或为空矩阵以表示应忽略其值。如果 n 大于计算 FFT 的维度,则 x 会被调整大小并用零填充。否则,如果 n 小于计算 FFT 的维度,则 x 会被截断。

如果使用三个参数调用,dim 是一个整数,指定执行 FFT 时所沿的矩阵维度。

另请参阅: ifft, fft2, fftn, fftw.

 
x = ifft (y)
x = ifft (y, n)
x = ifft (y, n, dim)

使用快速傅立叶变换(FFT)算法计算 y 的离散傅立叶逆变换。

逆 FFT 沿着数组的第一个非单一维度进行计算。因此,如果 y 是一个矩阵,ifft (y) 计算 y 每一列的逆 FFT。

如果使用两个参数调用,n 应为指定要使用的 y 元素数量的整数,或为空矩阵以表示应忽略其值。如果 n 大于计算逆 FFT 的维度,则 y 会被调整大小并用零填充。否则,如果 n 小于计算逆 FFT 的维度,则 y 会被截断。

如果使用三个参数调用,dim 是一个整数,指定执行逆 FFT 时所沿的矩阵维度。

另请参阅: fft, ifft2, ifftn, fftw.

 
B = fft2 (A)
B = fft2 (A, m, n)

使用快速傅立叶变换(FFT)算法计算 A 的二维离散傅立叶变换。

可选参数 mn 可用于指定要使用的 A 的行数和列数。如果其中任何一个大于 A 的大小,则 A 会被调整大小并用零填充。

如果 A 是一个多维矩阵,则 A 的每个二维子矩阵会被单独处理。

另请参阅: ifft2, fft, fftn, fftw.

 
A = ifft2 (B)
A = ifft2 (B, m, n)

使用快速傅立叶变换(FFT)算法计算 B 的二维离散傅立叶逆变换。

可选参数 mn 可用于指定要使用的 B 的行数和列数。如果其中任何一个大于 B 的大小,则 B 会被调整大小并用零填充。

如果 B 是一个多维矩阵,则 B 的每个二维子矩阵会被单独处理。

另请参阅: fft2, ifft, ifftn, fftw.

 
B = fftn (A)
B = fftn (A, size)

使用快速傅立叶变换(FFT)算法计算 A 的 N 维离散傅立叶变换。

可选参数 size 可用于指定要使用的 A 各维度的数量。如果 size 中任何维度的大小大于 A 中对应维度的大小,则 A 会被调整大小并用零填充。否则,A 会被截断。

另请参阅: ifftn, fft, fft2, fftw.

 
A = ifftn (B)
A = ifftn (B, size)

使用快速傅立叶变换(FFT)算法计算 B 的 N 维离散傅立叶逆变换。

可选参数 size 可用于指定要使用的 B 各维度的数量。如果 size 中任何维度的大小大于 B 中对应维度的大小,则 B 会被调整大小并用零填充。否则,B 会被截断。

另请参阅: fftn, ifft, ifft2, fftw.

 
method = fftw ("planner")
method = fftw ("planner", method)
wisdom = fftw ("wisdom")
wisdom = fftw ("wisdom", wisdom)
fftw ("dwisdom")

管理 FFTW 的计划程序(planner)设置和智慧(wisdom)。

从 FFTW 版本 3 开始,FFTW 库包含一种通过计划程序设置来加速 FFT 计算的方法。计划程序设置可以是以下之一:"estimate""measure""patient""exhaustive"。当通过 fftw ("planner", method) 设置方法时,所有后续的 FFT 计算都将使用该计划程序设置。默认值为 "estimate"

调用 fftw ("wisdom") 将以字符串形式返回当前的智慧。该字符串可以保存到磁盘并稍后使用 fftw ("wisdom", wisdom) 恢复。当智慧被恢复时,Octave 会强制将计划程序方法设置为 "measure",除非设置了 "estimate",以便复用存储的计划。FFTW 智慧在 Octave 会话之间不会自动保存。使用 fftw ("dwisdom") 可以将智慧转储到磁盘。

另请参阅: fft, ifft, fft2, ifft2, fftn, ifftn.

 
c = fftconv (x, y)
c = fftconv (x, y, n)

使用 FFT 对两个向量进行卷积。

c = fftconv (x, y) 返回长度为 length (x) + length (y) - 1 的向量。如果 xy 是两个多项式的系数向量,则返回值为乘积多项式的系数向量。

计算通过调用 fftfilt 函数使用 FFT。如果指定了可选参数 n,则使用 N 点 FFT。

另请参阅: deconv, conv, conv2.

 
y = fftfilt (b, x)
y = fftfilt (b, x, n)

使用 FFT 对 x 进行 FIR 滤波器 b 滤波。

如果 x 是一个矩阵,则对矩阵的每一列进行滤波。

如果提供了可选的第三个参数 nfftfilt 将使用重叠相加法(overlap-add method)通过 N 点 FFT 对 x 进行 b 滤波。FFT 大小必须是 2 的偶数幂,并且必须大于或等于 b 的长度。如果指定的 n 不满足这些条件,它将自动调整为满足条件的最近值。

另请参阅: filter, filter2.

 
y = filter (b, a, x)
[y, sf] = filter (b, a, x, si)
[y, sf] = filter (b, a, x, [], dim)
[y, sf] = filter (b, a, x, si, dim)

对数据 x 应用一维数字滤波器。

filter 返回以下线性时不变差分方程的解:

 N                   M
 SUM a(k+1) y(n-k) = SUM b(k+1) x(n-k)    for 1<=n<=length(x)
 k=0                 k=0

其中 N=length(a)-1,M=length(b)-1。 结果在 x 的第一个非单一维度上计算,如果提供了 dim 则在指定维度上计算。

等效形式的方程为:

          N                   M
 y(n) = - SUM c(k+1) y(n-k) + SUM d(k+1) x(n-k)  for 1<=n<=length(x)
          k=1                 k=0

其中 c = a/a(1),d = b/a(1)。

如果提供了第四个参数 si,它将被用作系统的初始状态,最终状态作为 sf 返回。状态向量是一个列向量,其长度等于最长系数向量的长度减一。如果未提供 si,初始状态向量将被设置为全零。

从 Z 变换的角度来看,y 是将离散时间信号 x 通过一个由以下有理系统函数表征的系统后的结果:

          M
          SUM d(k+1) z^(-k)
          k=0
H(z) = ---------------------
             N
        1 + SUM c(k+1) z^(-k)
            k=1

另请参阅: filter2, fftfilt, freqz.

 
y = filter2 (b, x)
y = filter2 (b, x, shape)

使用矩阵 b 对数据 x 进行二维 FIR 滤波。

如果可选参数 shape"full",则返回完整的卷积结果。默认值为 "same",返回与 x 相同大小的卷积中心部分。当 shape"valid" 时,仅返回未使用零填充边缘计算的卷积部分。

当使用 shape 的第三个参数时,根据 conv2 返回卷积的中心部分。这里 filter2 的作用与 conv2 (b, x) 相同,不过 b 相对于 x 有 180 度的旋转。

另请参阅: conv2, filter, fftfilt.

 
[h, w] = freqz (b, a, n)
[h, w] = freqz (b, a, n, "whole")
[h, f] = freqz (b, a, n, Fs)
[h, f] = freqz (b, a, n, Fs, "whole")
h = freqz (b, a, w)
[h, w] = freqz (b, a, n, "whole", Fs)

返回数字滤波器的复频率响应。

freqz 计算 Z 变换有理系统函数在其单位圆上的值。如果调用时不带输出参数,则会绘制幅度和相位响应。

有五种不同的调用方式:

freqz (b, a, n)

返回 [0, pi] 范围内 n 个点的复频率响应。

freqz (b, a, n, "whole")

返回 [0, 2*pi] 范围内 n 个点的复频率响应。

freqz (b, a, n, Fs)

返回 [0, Fs/2] 范围内 n 个点的复频率响应。

freqz (b, a, n, Fs, "whole")

返回 [0, Fs] 范围内 n 个点的复频率响应。

freqz (b, a, w)

返回在指定频率 w(弧度)下的复频率响应。

另请参阅: freqz_plot.

 
freqz_plot (w, h)
freqz_plot (w, h, freq_norm)

绘制由 freqz 返回的复频率响应。

如果可选参数 freq_norm 为 true(默认值),则频率向量 w 以归一化弧度为度量单位。如果 freq_norm 为 false 或未给出,则 w 以赫兹为度量单位。

另请参阅: freqz.

 
y = sinc (x)

计算 sinc 函数。

返回 sin (pi*x) / (pi*x)。

 
b = unwrap (x)
b = unwrap (x, tol)
b = unwrap (x, tol, dim)

通过适当地加减 2*pi 的倍数来展开弧度相位,以消除大于 tol 的跳变。

tol 默认为 pi。

unwrap 将沿着维度 dim 进行展开。如果未指定 dim,则默认为第一个非单一维度。

unwrap 忽略所有非有限输入值(Inf、NaN、NA)。

 
[a, b] = arch_fit (y, x, p, iter, gamma, a0, b0)

使用 Engle 原始 ARCH 论文中的评分算法,将 ARCH 回归模型拟合到时间序列 y

模型为:

y(t) = b(1) * x(t,1) + ... + b(k) * x(t,k) + e(t),
h(t) = a(1) + a(2) * e(t-1)^2 + ... + a(p+1) * e(t-p)^2

其中 e(t) 服从 N(0, h(t)),给定时间序列向量 y 到时间 t-1 和(普通)回归变量矩阵 xt。残差方差的回归阶数由 p 指定。

如果以 arch_fit (y, k, p) 的形式调用,其中 k 为正整数,则拟合 ARCH(k, p) 过程,即上述过程中 x 的第 t 行由下式给出:

[1, y(t-1), ..., y(t-k)]

可选地,可以指定评分算法的迭代次数 iter、更新因子 gamma 以及初始值 a0b0

 
y = arch_rnd (a, b, t)

模拟长度为 t 的 ARCH 序列,其中 AR 系数为 b,CH 系数为 a

结果 y(t) 遵循以下模型:

y(t) = b(1) + b(2) * y(t-1) + ... + b(lb) * y(t-lb+1) + e(t),

其中 e(t),给定 y 到时间 t-1,服从 N(0, h(t)),且

h(t) = a(1) + a(2) * e(t-1)^2 + ... + a(la) * e(t-la+1)^2
 
[pval, lm] = arch_test (y, x, p)

对于线性回归模型

y = x * b + e

执行拉格朗日乘数(LM)检验,检验无条件异方差的原假设对 CH(p) 的备择假设。

即模型为:

y(t) = b(1) * x(t,1) + ... + b(k) * x(t,k) + e(t),

给定 yt-1xte(t) 服从 N(0, h(t)),其中

h(t) = v + a(1) * e(t-1)^2 + ... + a(p) * e(t-p)^2,

原假设为 a(1) == ... == a(p) == 0。

如果第二个参数是标量整数 k,则在 k 阶线性自回归模型中进行相同的检验,即使用

[1, y(t-1), ..., y(t-k)]

作为 x 的第 t 行。

在原假设下,LM 近似服从自由度为 p 的卡方分布,pval 是检验的 p 值(1 - chi2cdf(lm, p))。

如果没有输出参数,pval 将被显示。

另请参阅: arch_fit, arch_rnd.

 
x = arma_rnd (a, b, v, t)

模拟长度为 t 的 ARMA 序列。

模型为:

x(n) = a(1) * x(n-1) + ... + a(la) * x(n-la) + e(n) + ...
       b(1) * e(n-1) + ... + b(lb) * e(n-lb)

其中 e(n) 是方差为 v 的白噪声。

 
m = autoreg_matrix (y, k)

给定一个时间序列(向量)y,返回一个矩阵,其中第 t 行包含 [1, y(t-1), ..., y(t-k)] 形式的自回归变量。

换句话说,对于 t = k+1, ..., T(其中 T = length(y)),m 的第 t 行为 [1, y(t-1), ..., y(t-k)]

结果的初始 k 行被设置为零。

 
c = bartlett (m)

返回长度为 m 的 Bartlett 窗口的滤波器系数。

Bartlett 窗口与三角窗口(triangular window)非常相似。Bartlett 窗口在两端始终为零,而三角窗口在两端非零。实际上,对于奇数长度 m,Bartlett 窗口是长度为 m+2 的三角窗口的中间 m 个点,且其两端被强制为零。

 
c = blackman (m)
c = blackman (m, "periodic")
c = blackman (m, "symmetric")

返回长度为 m 的 Blackman 窗口的滤波器系数。

如果给出了可选参数 "periodic",则返回窗口的周期形式。这相当于长度为 m+1 的窗口去掉最后一个系数。可选参数 "symmetric" 与不指定第二个参数等效。

关于 Blackman 窗口的定义,请参见例如 A.V. Oppenheim & R. W. Schafer, Discrete-Time Signal Processing

 
y = detrend (x, p)

如果 x 是一个向量,detrend (x, p) 从数据 x 中去除 p 阶多项式的最佳拟合。

如果 x 是一个矩阵,detrend (x, p)x 的每一列执行相同操作。

第二个参数 p 是可选的。如果未指定,则假定值为 1。这对应于去除线性趋势。

多项式的阶数也可以作为字符串给出,此时 p 必须是 "constant"(对应于 p=0)或 "linear"(对应于 p=1)。

另请参阅: polyfit.

 
[d, dd] = diffpara (x, a, b)

返回整合时间序列的差分参数估计值 d

使用 [2*pi*a/t, 2*pi*b/T] 范围内的频率进行估计。如果 b 省略,则使用 [2*pi/T, 2*pi*a/T] 范围。如果 ab 都省略,则使用 [2*pi/T, 2*pi*T/8] 范围。

输出 ddd 方差的一个估计值。

 
[la, v] = durbinlevinson (c, oldphi, oldv)

执行 Durbin-Levinson 算法的一步运算。

给定一个向量 c,其元素 [gamma_0, ..., gamma_t] 为自协方差,在给定 [gamma_0, ..., gamma_{t-1}][oldphi, oldv] 的情况下,返回用于预测 t-1 步的系数向量 la 和预测误差方差 v

如果 oldphioldv 省略,则执行完整的 Durbin-Levinson 算法。

 
y = fftshift (x)
y = fftshift (x, dim)

对向量 x 执行移位操作,用于 fftifft 函数,以便将频率 0 移动到向量或矩阵的中心。

如果 x 是一个包含 N 个元素的向量,对应 N 个以 dt 为间隔的时间样本,则 fftshift (fft (x)) 对应于频率

f = [ -(ceil((N-1)/2):-1:1), 0, (1:floor((N-1)/2)) ] * df

其中 df = 1 / (N * dt)

如果 x 是一个矩阵,对行和列同样成立。如果 x 是一个数组,则对每个维度同样成立。

可选参数 dim 可用于限制发生置换的维度。

另请参阅: ifftshift.

 
x = ifftshift (y)
x = ifftshift (y, dim)

撤销 fftshift 函数的操作。

对于偶数长度的 xfftshift 是其自身的逆操作,但奇数长度略有不同。

另请参阅: fftshift.

 
fd = fractdiff (x, d)

计算分数差分 (1-L)^d x,其中 L 表示滞后算子,d 大于 -1。

 
c = hamming (m)
c = hamming (m, "periodic")
c = hamming (m, "symmetric")

返回长度为 m 的 Hamming 窗口的滤波器系数。

如果给出了可选参数 "periodic",则返回窗口的周期形式。这相当于长度为 m+1 的窗口去掉最后一个系数。可选参数 "symmetric" 与不指定第二个参数等效。

关于 Hamming 窗口的定义,请参见例如 A.V. Oppenheim & R. W. Schafer, Discrete-Time Signal Processing

 
c = hanning (m)
c = hanning (m, "periodic")
c = hanning (m, "symmetric")

返回长度为 m 的 Hanning 窗口的滤波器系数。

如果给出了可选参数 "periodic",则返回窗口的周期形式。这相当于长度为 m+1 的窗口去掉最后一个系数。可选参数 "symmetric" 与不指定第二个参数等效。

关于 Hanning 窗口的定义,请参见例如 A.V. Oppenheim & R. W. Schafer, Discrete-Time Signal Processing

 
H = hurst (x)

估计时间序列 x 的 Hurst 参数。

 
pp = pchip (x, y)
yi = pchip (x, y, xi)

返回由点 xy 定义的分段三次 Hermite 插值多项式(Piecewise Cubic Hermite Interpolating Polynomial,PCHIP)。

如果使用两个参数调用,返回可用于 ppval 的分段多项式结构 pp,以在特定点处计算多项式值。

当使用第三个输入参数调用时,pchip 在点 xi 处计算 pchip 多项式的值。第三种调用形式等效于 ppval (pchip (x, y), xi)

变量 x 必须是长度为 n 的严格单调向量(递增或递减)。

y 可以是向量或数组。如果 y 是向量,则其长度必须与 x 相同 n。如果 y 是数组,则 y 的大小必须为 [s1, s2, …, sk, n]。该数组在内部被重塑为矩阵,其中前导维度由 s1 * s2 * … * sk 给出,然后该矩阵的每一行被单独处理。请注意,这与 interp1 正好相反,但这是为了与 MATLAB 兼容而做的。

另请参阅: spline, ppval, mkpp, unmkpp.

 
[Pxx, w] = periodogram (x)
[Pxx, w] = periodogram (x, win)
[Pxx, w] = periodogram (x, win, nfft)
[Pxx, f] = periodogram (x, win, nfft, Fs)
[Pxx, f] = periodogram (…, "range")
periodogram (…)

返回 x 的周期图(功率谱密度)。

可能的输入为:

x

数据向量。如果 x 是实值的,则估计单侧频谱。如果 x 是复值的,或者 "range" 指定为 "twosided",则估计全频谱。

win

窗口权重数据。如果窗口为空或未指定,则使用默认的矩形窗口。否则,在计算周期图之前,将窗口应用于信号(x .* win)。窗口数据必须是长度与 x 相同的向量。

nfft

频率箱的数量。默认值为 256 或大于 x 长度的下一个 2 的幂次(max (256, 2.^nextpow2 (length (x))))。如果 nfft 大于输入长度,则 x 会被补零至所需长度。如果省略 nfft,则将使用上述默认值。

Fs

采样频率。如果省略,则频率向量 w 以归一化弧度为度量单位。如果指定了 Fs,则频率向量 f 以赫兹为度量单位。

range

返回频谱的范围。"twosided" 返回全频谱,"onesided" 返回信号的一半频谱(适用于实值信号)。默认为 "onesided"

如果没有输出参数,将绘制周期图。

另请参阅: fft.

 
c = rectangular_sw (n, b)

矩形谱窗口。对于给定数量的数据 n 和带宽 b,返回窗口权重。

窗口名称用于搜索名为 win_sw 的函数。

 
y = sinetone (freq, rate, sec, ampl)

返回一个长度为 sec 秒的正弦波音调,采样率为 rate Hz,频率为 freq Hz,幅度为 ampl

参数 freqampl 是可选的,默认值为 440 Hz 和 1。

 
y = sinewave (m, n)

返回一个正弦波,其第一个 m 个点包含 n 个周期的正弦波。

 
sde = spectral_adf (c, win, b)

给定自协方差向量 c、窗口名称 win 和带宽 b,返回谱密度估计值。

窗口名称,例如 "triangle""rectangle",用于搜索名为 win_sw 的函数。

如果省略 win,则使用三角窗口。

如果省略 b,则使用 1 / sqrt (length (x))

另请参阅: spectral_xdf.

 
sde = spectral_xdf (x)
sde = spectral_xdf (x, win)
sde = spectral_xdf (x, win, b)

给定数据向量 x、窗口名称 win 和带宽 b,返回谱密度估计值。

窗口名称,例如 "triangle""rectangle",用于搜索名为 win_sw 的函数。

如果省略 win,则使用三角窗口。

如果省略 b,则使用 1 / sqrt (length (x))

另请参阅: spectral_adf.

 
savg = spencer (x)

返回 Spencer 的 15 点移动平均,应用于 x 的每一列。

 
y = stft (x)
y = stft (x, win_size)
y = stft (x, win_size, inc)
y = stft (x, win_size, inc, num_coef)
y = stft (x, win_size, inc, num_coef, win_type)
[y, c] = stft (…)

计算向量 x 的短时傅立叶变换,使用 num_coef 个系数,通过应用大小为 win_size 的数据窗口和 inc 个点的增量。

在计算傅立叶变换之前,将应用以下窗口之一:

"hanning"

win_type = 1

"hamming"

win_type = 2

"rectangle"

win_type = 3

窗口名称可以作为字符串传递,也可以通过 win_type 数字传递。

以下默认值用于未指定的参数:win_size = 80,inc = 24,num_coef = 64,win_type = 1。

y = stft (x, …) 根据 num_coef 个正频率返回傅立叶系数的绝对值。

[y, c] = stft (x, …) 返回整个 STFT 矩阵 y 和一个包含窗口大小、增量和窗口类型的 3 元素向量 c,这是 synthesis 函数所需要的。

另请参阅: synthesis.

 
x = synthesis (y, c)

根据信号的短时傅立叶变换 y 和指定窗口大小、增量和窗口类型的 3 元素向量 c 计算信号。

yc 的值可以通过以下方式得到:

[y, c] = stft (x , ...)

另请参阅: stft.

 
[a, v] = yulewalker (c)

给定自协方差向量 c [gamma_0, …, gamma_p],使用 Yule-Walker 估计拟合 AR(p) 模型。

返回 AR 系数 a 和白噪声方差 v