自学脚手架——“Data-Driven Science and Engineering” by steven L. brunton（二）

2.1 Fourier Series and Fourier Transforms

• 在P59中，公式（2.30）有错误应该为：

$$\hat{\text{f}}={\text{F}}_{2014}\text{f}=\left[\begin{array}{cc}{\text{I}}_{512}& {\text{D}}_{512}\\ {\text{I}}_{512}& -{\text{D}}_{512}\end{array}\right]\left[\begin{array}{cc}{\text{F}}_{512}& 0\\ 0& {\text{F}}_{512}\end{array}\right]\left[\begin{array}{c}{\text{f}}_{\mathrm{even}}\\ {\text{f}}_{\mathrm{odd}}\end{array}\right]$$

2.2 Discrete Fourier Transform (DFT) and Fast Fourier Transform (FFT)

• 在P63的Code 2.4中：

n = 128;
L = 30;
dx = L/(n);
x = -L/2:dx:L/2-dx; f = cos(x).*exp(-x.^2/25); % Function
df = -(sin(x).exp(-x.^2/25) + (2/25)x.f); % Derivative

%% Approximate derivative using finite Difference…
for kappa=1:length(df)-1
dfFD(kappa)=(f(kappa+1)-f(kappa))/dx;
end
dfFD(end+1) = dfFD(end);

%% Derivative using FFT (spectral derivative)
fhat = fft(f);
kappa = (2pi/L)[-n/2:n/2-1];
kappa = fftshift(kappa); % Re-order fft frequencies
dfhat = ikappa.*fhat;
dfFFT = real(ifft(dfhat));

%% Plotting commands
plot(x,df,’k’,’LineWidth’,1.5), hold on
plot(x,dfFD,’b–’,’LineWidth’,1.2)
plot(x,dfFFT,’r–’,’LineWidth’,1.2)
legend(’True Derivative’,’Finite Diff.’,’FFT Derivative’)

其中加粗的代码使用的是matlab中fftshift函数，下面截取了官方手册中的说明：

Y = fftshift(X) 通过将零频分量移动到数组中心，重新排列傅里叶变换 X。

• 如果 X 是向量，则 fftshift 会将 X 的左右两半部分进行交换。
• 如果 X 是矩阵，则 fftshift 会将 X 的第一象限与第三象限交换，将第二象限与第四象限交换。
• 如果 X 是多维数组，则 fftshift 会沿每个维度交换 X 的半空间。

Y = fftshift(X,dim) 沿 X 的维度 dim 执行运算。例如，如果 X 是矩阵，其行表示多个一维变换，则 fftshift(X,2) 会将 X 的每一行的左右两半部分进行交换。

平移向量元素
交换行向量的左右两半部分。如果一个向量的元素数为奇数，则中间的元素被视为属于向量的左半部分。

Xeven = [1 2 3 4 5 6];
fftshift(Xeven)

ans = 1×6

    4     5     6     1     2     3

Xodd = [1 2 3 4 5 6 7];
fftshift(Xodd)

ans = 1×7

    5     6     7     1     2     3     4

平移一维信号
分析信号的频率分量时，将零频分量平移到中心会很有帮助。
创建信号 S、计算其傅里叶变换，然后绘制功率。

fs = 100;               % sampling frequency
t = 0:(1/fs):(10-1/fs); % time vector
S = cos(2*pi*15*t);
n = length(S);
X = fft(S);
f = (0:n-1)*(fs/n);     %frequency range
power = abs(X).^2/n;    %power
plot(f,power)

平移零频分量，然后绘制以零为中心的功率。

Y = fftshift(X);
fshift = (-n/2:n/2-1)*(fs/n); % zero-centered frequency >range
powershift = abs(Y).^2/n;     % zero-centered power
plot(fshift,powershift)

对应的还有[ifftshift函数]

2.4 Gabor Transform and the Spectrogram

• 在Code 2.13 中：

t = 0:0.001:2;
f0 = 50;
f1 = 250;
t1 = 2;
x = chirp(t,f0,t1,f1,’quadratic’);
x = cos(2pit.*(f0 + (f1-f0)*t.^2/(3*t12)));
% There is a typo in Matlab documentation…
% … divide by 3 so derivative amplitude matches frequency
spectrogram(x,128,120,128,1e3,’yaxis’)

其中加粗的代码使用了matlab中的spectrogram函数，下面截取了官方的使用说明：

• 在Figure 2.21 中：

其中左侧的图应该顺时针旋转，即与Figure 2.20中的功率谱密度一致。

• 在Code 2.14中

% Download mp3read from http://www.mathworks.com/matlabcentral/
fileexchange/13852-mp3read-and-mp3write
[Y,FS,NBITS,OPTS] = mp3read(’beethoven.mp3’);

%% Spectrogram using ‘spectrogram’ comand
T = 40; % 40 seconds
y=Y(1:T*FS); % First 40 seconds
spectrogram(y,5000,400,24000,24000,’yaxis’);

%% Spectrogram using short-time Fourier transform ‘stft’
wlen = 5000; % Window length
h=400; % Overlap is wlen - h
% Perform time-frequency analysis
[S,f,t_stft] = stft(y, wlen, h, FS/4, FS); % y axis 0-4000HZ

imagesc(log10(abs(S))); % Plot spectrogram (log-scaled)

其中使用了matlab中的stft函数，下面截取官方文档中的说明：

• Mexican hat wavelet

Mexican hat wavelet也被称为Ricker wavelet，公式为：

$$\psi (t)=(1-t^2)e^{-t^2/2}$$

在matlab中有现成的函数 mexihat表示：

lb = -5;
ub = 5;
N = 1000;
[psi,xval] = mexihat(lb,ub,N);
plot(xval,psi)
title('Mexican Hat Wavelet')

如图所见，几乎100%的能量集中在[-5,5]的区间中，如果需要更高精度，一般可以取[-8,8]。

注意高斯概率密度函数的二阶导数与这个公式成正比关系（差一个系数）。

2.6 2D Transforms and Image Processing

为何在Code 2.15中做每行的傅里叶变换的时候没有用spectrogram这个函数呢？fft2是用的哪一种？