本章介绍 Octave 中可用的信号处理和快速傅立叶变换函数。快速傅立叶变换的计算使用 FFTW 或 FFTPACK 库,具体取决于 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 时所沿的矩阵维度。
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 时所沿的矩阵维度。
B = fft2 (A) ¶B = fft2 (A, m, n) ¶使用快速傅立叶变换(FFT)算法计算 A 的二维离散傅立叶变换。
可选参数 m 和 n 可用于指定要使用的 A 的行数和列数。如果其中任何一个大于 A 的大小,则 A 会被调整大小并用零填充。
如果 A 是一个多维矩阵,则 A 的每个二维子矩阵会被单独处理。
A = ifft2 (B) ¶A = ifft2 (B, m, n) ¶使用快速傅立叶变换(FFT)算法计算 B 的二维离散傅立叶逆变换。
可选参数 m 和 n 可用于指定要使用的 B 的行数和列数。如果其中任何一个大于 B 的大小,则 B 会被调整大小并用零填充。
如果 B 是一个多维矩阵,则 B 的每个二维子矩阵会被单独处理。
B = fftn (A) ¶B = fftn (A, size) ¶使用快速傅立叶变换(FFT)算法计算 A 的 N 维离散傅立叶变换。
可选参数 size 可用于指定要使用的 A 各维度的数量。如果 size 中任何维度的大小大于 A 中对应维度的大小,则 A 会被调整大小并用零填充。否则,A 会被截断。
A = ifftn (B) ¶A = ifftn (B, size) ¶使用快速傅立叶变换(FFT)算法计算 B 的 N 维离散傅立叶逆变换。
可选参数 size 可用于指定要使用的 B 各维度的数量。如果 size 中任何维度的大小大于 B 中对应维度的大小,则 B 会被调整大小并用零填充。否则,B 会被截断。
method = fftw ("planner") ¶method = fftw ("planner", method) ¶wisdom = fftw ("wisdom") ¶wisdom = fftw ("wisdom", wisdom) ¶("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") 可以将智慧转储到磁盘。
c = fftconv (x, y) ¶c = fftconv (x, y, n) ¶使用 FFT 对两个向量进行卷积。
c = fftconv (x, y) 返回长度为 length (x) + length (y) - 1 的向量。如果 x 和 y 是两个多项式的系数向量,则返回值为乘积多项式的系数向量。
计算通过调用 fftfilt 函数使用 FFT。如果指定了可选参数 n,则使用 N 点 FFT。
y = fftfilt (b, x) ¶y = fftfilt (b, x, n) ¶使用 FFT 对 x 进行 FIR 滤波器 b 滤波。
如果 x 是一个矩阵,则对矩阵的每一列进行滤波。
如果提供了可选的第三个参数 n,fftfilt 将使用重叠相加法(overlap-add method)通过 N 点 FFT 对 x 进行 b 滤波。FFT 大小必须是 2 的偶数幂,并且必须大于或等于 b 的长度。如果指定的 n 不满足这些条件,它将自动调整为满足条件的最近值。
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
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 度的旋转。
[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.
(w, h) ¶(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 和(普通)回归变量矩阵 x 到 t。残差方差的回归阶数由 p 指定。
如果以 arch_fit (y, k, p) 的形式调用,其中 k 为正整数,则拟合 ARCH(k, p) 过程,即上述过程中 x 的第 t 行由下式给出:
[1, y(t-1), ..., y(t-k)]
可选地,可以指定评分算法的迭代次数 iter、更新因子 gamma 以及初始值 a0 和 b0。
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),
给定 y 到 t-1 和 x 到 t,e(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 将被显示。
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] 范围。如果 a 和 b 都省略,则使用 [2*pi/T, 2*pi*T/8] 范围。
输出 dd 是 d 方差的一个估计值。
[la, v] = durbinlevinson (c, oldphi, oldv) ¶执行 Durbin-Levinson 算法的一步运算。
给定一个向量 c,其元素 [gamma_0, ..., gamma_t] 为自协方差,在给定 [gamma_0, ..., gamma_{t-1}] 和 [oldphi, oldv] 的情况下,返回用于预测 t-1 步的系数向量 la 和预测误差方差 v。
如果 oldphi 和 oldv 省略,则执行完整的 Durbin-Levinson 算法。
y = fftshift (x) ¶y = fftshift (x, dim) ¶对向量 x 执行移位操作,用于 fft 和 ifft 函数,以便将频率 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 函数的操作。
对于偶数长度的 x,fftshift 是其自身的逆操作,但奇数长度略有不同。
另请参阅: 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) ¶返回由点 x 和 y 定义的分段三次 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 兼容而做的。
[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") ¶(…) ¶返回 x 的周期图(功率谱密度)。
可能的输入为:
数据向量。如果 x 是实值的,则估计单侧频谱。如果 x 是复值的,或者 "range" 指定为 "twosided",则估计全频谱。
窗口权重数据。如果窗口为空或未指定,则使用默认的矩形窗口。否则,在计算周期图之前,将窗口应用于信号(x .* win)。窗口数据必须是长度与 x 相同的向量。
频率箱的数量。默认值为 256 或大于 x 长度的下一个 2 的幂次(max (256, 2.^nextpow2 (length (x))))。如果 nfft 大于输入长度,则 x 会被补零至所需长度。如果省略 nfft,则将使用上述默认值。
采样频率。如果省略,则频率向量 w 以归一化弧度为度量单位。如果指定了 Fs,则频率向量 f 以赫兹为度量单位。
返回频谱的范围。"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。
参数 freq 和 ampl 是可选的,默认值为 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 计算信号。
y 和 c 的值可以通过以下方式得到:
[y, c] = stft (x , ...)
另请参阅: stft.
[a, v] = yulewalker (c) ¶给定自协方差向量 c [gamma_0, …, gamma_p],使用 Yule-Walker 估计拟合 AR(p) 模型。
返回 AR 系数 a 和白噪声方差 v。