傅里叶变换及其应用

前言

最近写智能设计与建造的作业，偶然兴起想把地震动响应(时域信号)转换成频域信号，需要用到傅里叶变换，于是就去系统学了学，做了相关笔记(手稿)，后面有兴趣还会做一些简单的项目实战一下。
参考资料：

  1. 李永乐讲傅里叶变换: 只推荐前半部分，后半部分就在乱讲。
  2. 傅里叶变换通俗解释: 对于三角函数系的正交表示提了一嘴，有利于我们记忆。
  3. 从傅里叶级数到傅里叶变换：详细的数学推导
  4. 傅里叶变换推导: 这个视频很好地弥补了李永乐老师讲的不足的地方。
  5. 本篇博客封面：本篇博客封面引用了该文章里的图片，内容并未详细阅读。

理论推导

实战

本文选取的地震波文件:acc.txt，使用离散傅里叶变换将地震加速度时程曲线转化成频域曲线，相关代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
clear all
clc


file_acc = fopen('./acc.txt','r');
acc_data = fscanf(file_acc,'%f',[1,inf]);
acc_data = reshape(acc_data,[],1);
acc_data = 9.8 * acc_data;  % Convert to m/s^2
fclose(file_acc);

rows_acc_data = size(acc_data,1);
total_time = (rows_acc_data-1) * 0.01;  % Sample interval is 0.01s
time_data = 0:0.01:total_time;
time_acc = [time_data', acc_data];

%% Manual DFT
base_freq = 2 * pi / (rows_acc_data-1);
Freq = zeros(rows_acc_data-1, 2);
for k = 1:rows_acc_data-1
    Freq(k,1) = (k-1) / total_time;
    temp_freq = 0;
    for n = 1:rows_acc_data-1
        temp_freq = temp_freq + time_acc(n,2) * exp(-1i * (n-1) * (k-1) * 2 * pi / (rows_acc_data-1));
    end
    Freq(k,2) = abs(temp_freq);
end

n_f = floor((rows_acc_data-1)/2);  % Nyquist frequency index

%% MATLAB Built-in FFT
Fs = 100;  % Sampling frequency (1/0.01)
L = rows_acc_data-1;
Y = fft(acc_data(1:end-1));
P2 = abs(Y);
P1 = P2(1:floor(L/2)+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs * (0:(L/2))/L;

%% Plot comparison
figure(1);
% Manual DFT plot
subplot(2,1,1);
plot(Freq(1:n_f,1), Freq(1:n_f,2), 'b-', 'LineWidth', 1);
title("Manual DFT");
xlabel("Frequency (Hz)");
ylabel("Magnitude");
grid on;

% MATLAB FFT plot
subplot(2,1,2);
plot(f, P1, 'r-', 'LineWidth', 1);
title("MATLAB FFT");
xlabel("Frequency (Hz)");
ylabel("Magnitude");
grid on;

最后结果也是和matlab中fft计算结果吻合较好。
加速度时程曲线：
加速度时程曲线

利用DFT公式计算结果与MATLAB中FFT计算结果对比：
DFT与FFT