傅里叶变换

傅立叶变换是用于将一个数字向量转换为另一个数字向量的公式。第一向量由在个别时间点取的一组信号值组成。第二矢量由表征谐波信号的幅度和频率的离散值集合组成，从中可以通过在时间或空间上重构所观察到的信号来构成。

考虑向量，其中包括以相同的时间间隔进行测量。这样的矢量的傅立叶变换描述如下:

角频率是一个复杂的根единицы, -假想单位。对于两个向量和 ,索引和它们取值从0到 .

功能 fft 从图书馆 FFTW 使用快速傅立叶变换（FFT）算法将矢量分解为傅立叶级数。考虑正弦信号，各自以相等的时间间隔取下。测量的时刻形成一个向量 . 该信号由频率为15和20Hz的两个分量组成。让我们研究信号，每1/50秒测量10秒。

Pkg.add(["WAV"])

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

将信号分解为傅立叶级数并创建矢量，其离散分量在频域中形成信号的样本。

using FFTW

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

让我们绘制变换的结果，其中谐波分量的频率从横坐标导出，幅度从纵坐标导出。

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

该图显示了四个脉冲串，尽管原始信号有两个频率为15和20Hz的分量。事实是，图表的右半部分反映了左半部分。信号在时域中的离散傅立叶变换也具有其自身的周期性–矢量的前半部分对应于"正频率值"的区域，后半部分对应于分量的"负频率值"的区域。 30和35Hz区域中的突发实际上对应于值-20和-15Hz。

在实践中，通常只使用傅立叶级数的一半而不执行不必要的计算就足够了。为此，您可以使用该函数 rfft. 但为了清楚起见，我们将输出整个傅立叶级数，使用函数将右半部分移位到负频率区域 fftshift.

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

有噪声信号

在现实场景中，测量信号通常包括有用分量和随机噪声，这使得无法清楚地区分谐波分量。使用傅立叶变换可以帮助切断随机分量，只留下系统的振荡过程。考虑一个新的测量向量 xnoise，信号相加的结果 x 和随机高斯噪声。

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

在实践中，信号功率对频率的依赖比振幅对频率的依赖更经常地使用。功率等于傅立叶级数的分量的幅度除以级数中的点数的平方。让我们计算并绘制信号功率谱，轴中心为零频率 . 尽管有噪声，但我们可以清楚地看到信号功率图上的两个尖峰。

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

转换率

要计算整个系列的傅立叶变换从向量元素使用初始公式，您需要执行以下操作用浮点数计算。快速傅立叶变换算法需要以下数量级这使得这种方法非常有吸引力，因为节省了时间，特别是如果正在处理的信号由几百万个点组成。而且，FFT算法的特定实现可以在若干特定情况下甚至更快地计算。例如，在当在某种程度上等于2。

考虑一个音频文件 bluewhale.wav -在水下麦克风上录制蓝鲸的声学通信。

using WAV
x, fs = wavread( "bluewhale.wav" );

蓝鲸在次声波中进行通信，声学振动的频带几乎不与人类听觉敏感的频率重叠。因此，所呈现的数据中的时间矢量被压缩10的因子，以便音频信号的音高对我们的感知是可访问的。

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

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 );

可视化此信号:

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

现在让我们将信号长度增加到最接近的数字。在上升方向上。我们会得到一段较长的信号。并为信号的这一新段计算傅立叶级数。如果信号长度不是从系列中提取的（当），那么函数 fft 最有可能的是，无论如何，我都会用零补充信号到最近的长度。 . 这样做是为了显着加快FFT算法，特别是如果点数的值没有小的除数。

m = length( whaleMoan ); 
n = nextpow(2, m);

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

y = fft( whaleMoan_zeros );

让我们绘制这个信号的功率谱。根据图表，声学信号包括频率约为17Hz的主要谐波-声音以及其他几个明显的谐波分量。

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 )

正弦曲线的相位检测

使用傅立叶变换，也可以构造所研究信号的相位谱。为了举例，我们将创建一个信号，其中频率为15和20Hz的两个谐波信号组合在一起。第一个信号将是具有相位值的余弦，并且第二信号将是没有相移的正弦波。信号的采样频率为100Hz，信号长度为1s。

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

将信号分解为傅立叶级数，并将幅度绘制为频率的函数。

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 )

现在让我们看看傅立叶级数的系数的相位部分，同时去除具有低幅度的值。该图显示了信号的相位谱作为频率的函数。

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 )

结论

我们已经构建了几个不同的频谱（振幅，功率和相位）的基础上分解的信号到傅立叶级数。该功能帮助我们做到了这一点。 fft，使用FFT算法，以及其他几个函数（例如 fftshift因此，我们能够直观地显示我们感兴趣的输入信号的所有属性。