Engee 文档
Notebook

傅立叶变换

傅立叶变换*是一个用于将一个数字向量转换成另一个数字向量的公式。第一个矢量由一组在离散时间时刻提取的信号值组成。第二个矢量由一组离散的数值组成，这些数值表征谐波信号的振幅和频率，通过在时间或空间上对其进行重构，可以得到观测到的信号。

考虑一个矢量$x$ ，其中包括以相等时间间隔进行的$n$ 测量。该矢量的傅立叶变换描述如下：

$$y_{k+1} = \sum_{j=0}^{n-1}\omega^{jk}x_{j+1}$$

角频率$\omega = e^{-2\pi i/n}$ 是复数1 的根，$i$ 是虚数单位。对于$x$ 和$y$ 这两个向量，索引$j$ 和$k$ 的取值范围为 0 至$n-1$ 。

函数库FFTW 中的函数fft 使用快速傅里叶变换 (FFT) 算法将矢量分解为傅里叶级数。考虑正弦信号$x$ ，每次测量的时间间隔相等。测量时刻构成一个矢量$t$ 。信号由频率为 15 和 20 Hz 的两个分量组成。让我们研究一下在 10 秒内每隔 1/50 秒测量一次的信号。

In [ ]: 
Pkg.add(["WAV"])
In [ ]: 
Ts = 1/50;
t = 0:Ts:10-Ts;                     
x = sin.(2pi * 15t) .+ sin.( 2pi * 20t );
plot( t, x, xlabel="Время (с)", ylabel="Амплитуда", leg=false )
Out[0]:

将信号分解为傅里叶级数，并创建一个矢量$f$ ，其离散分量构成频域中的信号样本。

In [ ]: 
using FFTW

y = fft(x);
fs = 1/Ts;
f = (0:length(y)-1)*fs/length(y);

让我们绘制一下变换结果，其中的横座标表示谐波分量的频率，纵座标表示振幅。

In [ ]: 
plot( f, abs.(y), title="Амплитудный спектр", xlabel="Частота (Hz)", ylabel="Амплитуда" )
Out[0]:

虽然原始信号有两个频率分别为 15 赫兹和 20 赫兹的分量，但图中显示的是四个脉冲串。重点是图表的右半部分反映了左半部分。信号在时域中的离散傅里叶变换也有其周期性--矢量的前半部分对应于 "正频率值 "区域，后半部分对应于各分量的 "负频率值 "区域。30 和 35 Hz 区域的脉冲串实际上对应于 -20 和 -15 Hz 的值。

在实际应用中，通常只需使用傅立叶级数的一半即可，无需进行不必要的计算。为此，您可以使用函数rfft 。不过，为了清楚起见，我们还是通过使用函数fftshift 将右半部分移至负频率区域来推导整个傅里叶级数。

In [ ]: 
n = length(x);
fshift = (-n/2:n/2-1)*(fs/n);
yshift = fftshift(y);
plot( fshift, abs.(yshift), title="Амплитудный спектр", xlabel="Частота (Uw)", ylabel="Амплитуда" )
Out[0]:

噪声信号

在现实世界中，测量信号往往包含有用成分和随机噪声，因此无法明确区分谐波成分。使用傅立叶变换可以帮助去除随机成分，只留下系统的振荡过程。考虑一个新的测量矢量xnoise ，它是信号x 和随机高斯噪声相加的结果。

In [ ]: 
using Random
Random.seed!(0)
xnoise = x .+ 2.5 .* randn( size(t) );

在实际应用中，信号功率的频率相关性比振幅的频率相关性更常用。功率等于傅里叶级数分量的振幅平方除以级数中的点数。让我们计算并绘制信号的功率谱，零频率位于轴中心$x$ 。尽管存在噪声，我们仍能在信号功率图上清楚地看到两个脉冲串。

In [ ]: 
ynoise = fft( xnoise );
ynoiseshift = fftshift( ynoise );    
power = abs.( ynoiseshift ).^2 ./ n; 
plot( fshift, power, title="Спектр мощности", xlabel="Частота (Гц)", ylabel="Мощность сигнала" )
Out[0]:

转换率

使用原始公式计算$n$ 向量$y$ 的整个系列元素的傅里叶变换，大约需要$n^2$ 浮点运算。快速傅立叶变换算法大约需要$n log n$ 次运算，这使得这种方法非常有吸引力，因为可以节省时间，尤其是当要处理的信号由几百万个点组成时。此外，在一些特定情况下，FFT 算法的特定实现可以更快地计算。例如，在$n$ 在某种程度上是 2 的情况下。

考虑音频文件bluewhale.wav --蓝鲸通过水下麦克风录制的声学通讯记录。

In [ ]: 
using WAV
x, fs = wavread( "bluewhale.wav" );

蓝鲸用次声波进行交流，声波振荡的频段与人类听觉敏感的频率几乎不重叠。因此，所提供数据中的时间矢量被压缩了 10 倍，以便我们能够感知声音信号的高度。

下一个单元允许您回放从文件中读取的声音信号。

In [ ]: 
using Base64

buf = IOBuffer();
wavwrite(x, buf; Fs=fs);
data = base64encode(unsafe_string(pointer(buf.data), buf.size));
markup = """<audio controls="controls" {autoplay}>
              <source src="data:audio/wav;base64,$data" type="audio/wav" />
              Your browser does not support the audio element.
              </audio>"""
display( "text/html", markup );

让我们将信号可视化：

In [ ]: 
whaleMoan = x[24500:31000];
t = 10 * (0:1/fs:(length(whaleMoan))/fs);
plot( t, whaleMoan, xlabel="Время (с)", ylabel="Амплитуда", label=false, xlimits=(0,t[end]) )
Out[0]:

现在，让我们将信号的长度向上增加到最接近的数字$2^N$ 。我们将得到一段更长的信号。然后计算这个新信号段的傅立叶级数。如果信号长度没有取自序列$2^N$ （在$N = 1,2\dots$ 处），函数fft ，很可能仍会在信号中增加零，达到最接近的长度$2^N$ 。这样做是为了大大加快 FFT 算法的速度，尤其是在点数值没有小除数的情况下。

In [ ]: 
m = length( whaleMoan ); 
n = nextpow(2, m);

whaleMoan_zeros = zeros( n )
whaleMoan_zeros[1:length(whaleMoan)] = whaleMoan

y = fft( whaleMoan_zeros );

让我们绘制该信号的功率谱图。从图中可以看出，声音信号包括主谐波--频率约为 17 赫兹的声音，以及其他几个明显的谐波成分。

In [ ]: 
f = (0:n-1) .* (fs/n)/10; # вектор значений частоты
power = abs.(y).^2/n;     # спектр мощности
plot( f[1:Int(floor(n/2))],power[1:Int(floor(n/2))],
    xlabel = "Частота (Гц)", ylabel = "Мощность", leg=false )
Out[0]:

确定正弦波的相位

傅立叶变换也可用于构建所研究信号的相位频谱。在此示例中，我们将创建一个信号，将频率为 15 和 20 Hz 的两个谐波信号相加。第一个信号是相位值为$-\pi/4$ 的余弦波，第二个信号是没有相位偏移的正弦波。信号的采样频率为 100 赫兹，信号长度为 1 秒。

In [ ]: 
fs = 100;
t = 0:1/fs:1-1/fs;
x = @. cos(2*pi*15*t - pi/4) - sin(2*pi*40*t);

让我们对信号进行傅里叶级数分解，并绘制振幅与频率的函数关系图。

In [ ]: 
y = fft(x);
z = fftshift(y);

ly = length(y);
f = (-ly/2:ly/2-1)/ly*fs;

plot( f, abs.(z), st=:stem, marker=(5, :white, :o),
    markerstrokecolor=1,
    xlabel="Частота (Гц)", ylabel="|y|", grid=true )
Out[0]:

现在我们来看看傅里叶级数系数的相位部分，去除振幅较低的值。该图显示了信号的相位频谱与频率的函数关系。

In [ ]: 
tol = 1e-6;
z[abs.(z) .< tol] .= 0;

theta = angle.(z);

plot( f, theta/pi, st=:stem, marker=(5, :white, :o),
    markerstrokecolor=1,
    xlabel="Частота (Гц)", ylabel = "Фаза / π", grid=true )
Out[0]:

结论

我们根据信号的傅里叶级数展开构建了多个不同的频谱（振幅、功率和相位）。使用 FFT 算法的函数fft 和其他几个函数（如fftshift ）帮助我们完成了这项工作。因此，我们能够清楚地显示我们感兴趣的输入信号的所有属性。