Matlab Math
Cleve Morler著
陈文斌（ wbchen@fudan.edu.cn)
复旦大学 2002
Fourier 变换
Touch-tone Dialing
我们每天都在使用 Fourier变换，比如移动电话,CD,DVD,JPEG
1
ABC
2
DEF
3
GHI
4
JKL
5
MNO
6
PRS
7
TUV
8
WXY
9
*
OPER
0 #
697
770
852
941
1209 1336 1477
Dtmfanalyze.m
Dtmfsynthesize.m
Finite Fourier Transform
??
?
?? ?
1
0
11
n
j
j
jk
k yY ?
nie /2?? ??
离散 Fourier变换
FyY ? jkjkf ???? 1,1
nIFF H ?
YFny H1? ??
?
?? ?
1
0
11
1 n
j
jk
kj Yny ?
nie /2?? ??
离散反 Fourier变换
omega=exp(-2*pi*i/n); j=0:n-1;
k=j'; F=omega.^(k*j);
F=fft(eye(n));
fftgui
fftgui
fftgui
fftgui
Nyquist point
fftgui
对称性
?? jyY )(r e a l 1
如果 y是长度为 n的实向量,Y=fft(y)
0)(i m a g 1 ?Y
12/,...,0),(r e a l)(r e a l 2 ??? ?? njYY jnj
12/,...,0),(i m a g)(i m a g 2 ???? ?? njYY jnj
Sunspots
For centuries people have noted that the face of the sun is not
constant or uniform in appearance,but that correlated with weather
and other economically significant terrestrial phenomena,In 1848,
Rudolf Wolfer proposed a rule that combined the number and size
of these sunspots into a single index,Using archival records,
astronomers have applied Wolfer's rule to determine sunspot
activity back to the year 1700,Today the sunspot index is measured
by many astronomers and the worldwide distributioin of the data is
coordinated by the Sunspot Index Data Center at Royal
Observatory of Belgium.
http://sidc.oma.be/html/sunspot.html
http://www.noao.edu/image_gallery/html/im0404.html
load sunspot.dat; t = sunspot(:,1)';
wolfer = sunspot(:,2)'; n=length(wolfer);
Sunspots
c=polyfit(t,wolfer,1); trend=polyval(c,t);
plot(t,[wolfer;trend],'-',t,wolfer,'k.');
xlabel('year');ylabel('Wolfer index');
title('Sunspot index with linear trend');
1700 1750 1800 1850 1900 1950 20000
20
40
60
80
100
120
140
160
180
200
year
W
olf
er
ind
ex
最高和最低好像蛮有规律的呀
试试 FFT吧！
y=wolfer-trend; Y=fft(y); Y(1)=[];
pow=abs(Y(1:n/2)).^2; pmax=20e6;
f=(1:n/2)/n;
plot([f; f],[0*pow; pow],'c-',f,pow,'b.',…
'linewidth',2,'markersize',16);
axis([0,5 0 pmax]); xlabel('cycles/year');
ylabel('power'); title('Periodogram');
Sunspots
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.50
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2 x 10
7
cycles/year
po
we
r
Periodogram
好高噢！
k=1:36;pow=pow(k);
ypk=n./k(2:2:end); %Years per cycle
plot([k;k],[0*pow;pow],'c-',k,pow,'b.',…
'linewidth',2,'markersize',16)
axis([0 max(k)+1 0 pmax]);
set(gca,'xtick',k(2:2:end));
xticklabels=sprintf('%5.1f|',ypk);
set(gca,'xticklabel',xticklabels);
xlabel('years/cycle');ylabel('power')
title('Periodogram')
Sunspots
144.0 72.0 48.0 36.0 28.8 24.0 20.6 18.0 16.0 14.4 13.1 12.0 11.1 10.3 9.6 9.0 8.5 8.00
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2 x 10
7
years/cycle
po
we
r
Periodogram
原来如此
Fast Finite Fourier Transform
??
?
?? ???
1
0
11 1,.,,,0,
n
j
j
jk
k nkyY ?
每个 k需要 n个乘 n个加，总共 2n× n次浮点运算
（不包括矩阵的生成）
给出你计算机的频率，算一下需要多长时间？
当代的 FFT,1965
Jonh Tukey(Princeton University)
John Cooley(IBM Research)
)()l o g( 22 nOnnO ?
nn 2lo g3
Fastest Fourier Transform in
the West
试一下，你的机器运行 1024*1024的 fft运算要多长时间
y=rand(1024); Y=fft(y); plot(y,'.');
figure; plot(Y,'.');
Fastest Fourier Transform in the West
http://www.fftw.org MIT
Fast Finite Fourier Transform
???? ? s i nc o s/2 ie nin ???? ?
nn ?? ?22
??
???
?
?
?
?
?
?
??
?
?
??
??
???
12/
0
22
2
12/
0
12
2
odd
1
e v e n
1
1
0
11
n
j
j
jkk
n
j
j
jk
j
j
jk
j
j
jk
n
j
j
jk
k
yy
yyyY
???
???
如果 n是偶数且 k≤n/2-1,把整个求和分成两个部分
Fast Finite Fourier Transform
n=length(y);
omega=exp(-2*pi*i/n); k=(0:n/2-1)';
w=omega.^k; u=fft(y(1:2:n-1));
v=w.*fft(y(2:2:n));
fft(y)=[u+v; u-v];
如果 n是 2的 p次，这个过程可以重复 p次，每次的计算
量是 O(n),总的计算量是
)l o g()( 2 nnOnpO ?
ffttx.m
Fast Finite Fourier Transform
ffttx.m
原型：
fftmatrixgui.m
图的稀疏性，影响 FFT算法的速度
F=fft(eye(n,n)); plot(fft(eye(n,n)))
Fourier Transforms and Series
??
?
?? ???
1
0
11 1,.,,,0,
n
j
j
jk
k nkyY ?
1,...,0,1
1
0
11 ??? ?
?
?
?? njYny
n
j
jk
kj ?
? ??? ?? dtetfF ti??? 2)()(
? ???? ?? ?? deFtf ti2)()(
如果 f是周期函数,f(t+L)=f(t)
?? ? ?
?
???
?? 2/
2/
2/2 )(1,)( L
L
i j t
j
j
Li j t
j dtetfLcectf
??
Cleve Morler著
陈文斌（ wbchen@fudan.edu.cn)
复旦大学 2002
Fourier 变换
Touch-tone Dialing
我们每天都在使用 Fourier变换，比如移动电话,CD,DVD,JPEG
1
ABC
2
DEF
3
GHI
4
JKL
5
MNO
6
PRS
7
TUV
8
WXY
9
*
OPER
0 #
697
770
852
941
1209 1336 1477
Dtmfanalyze.m
Dtmfsynthesize.m
Finite Fourier Transform
??
?
?? ?
1
0
11
n
j
j
jk
k yY ?
nie /2?? ??
离散 Fourier变换
FyY ? jkjkf ???? 1,1
nIFF H ?
YFny H1? ??
?
?? ?
1
0
11
1 n
j
jk
kj Yny ?
nie /2?? ??
离散反 Fourier变换
omega=exp(-2*pi*i/n); j=0:n-1;
k=j'; F=omega.^(k*j);
F=fft(eye(n));
fftgui
fftgui
fftgui
fftgui
Nyquist point
fftgui
对称性
?? jyY )(r e a l 1
如果 y是长度为 n的实向量,Y=fft(y)
0)(i m a g 1 ?Y
12/,...,0),(r e a l)(r e a l 2 ??? ?? njYY jnj
12/,...,0),(i m a g)(i m a g 2 ???? ?? njYY jnj
Sunspots
For centuries people have noted that the face of the sun is not
constant or uniform in appearance,but that correlated with weather
and other economically significant terrestrial phenomena,In 1848,
Rudolf Wolfer proposed a rule that combined the number and size
of these sunspots into a single index,Using archival records,
astronomers have applied Wolfer's rule to determine sunspot
activity back to the year 1700,Today the sunspot index is measured
by many astronomers and the worldwide distributioin of the data is
coordinated by the Sunspot Index Data Center at Royal
Observatory of Belgium.
http://sidc.oma.be/html/sunspot.html
http://www.noao.edu/image_gallery/html/im0404.html
load sunspot.dat; t = sunspot(:,1)';
wolfer = sunspot(:,2)'; n=length(wolfer);
Sunspots
c=polyfit(t,wolfer,1); trend=polyval(c,t);
plot(t,[wolfer;trend],'-',t,wolfer,'k.');
xlabel('year');ylabel('Wolfer index');
title('Sunspot index with linear trend');
1700 1750 1800 1850 1900 1950 20000
20
40
60
80
100
120
140
160
180
200
year
W
olf
er
ind
ex
最高和最低好像蛮有规律的呀
试试 FFT吧！
y=wolfer-trend; Y=fft(y); Y(1)=[];
pow=abs(Y(1:n/2)).^2; pmax=20e6;
f=(1:n/2)/n;
plot([f; f],[0*pow; pow],'c-',f,pow,'b.',…
'linewidth',2,'markersize',16);
axis([0,5 0 pmax]); xlabel('cycles/year');
ylabel('power'); title('Periodogram');
Sunspots
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.50
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2 x 10
7
cycles/year
po
we
r
Periodogram
好高噢！
k=1:36;pow=pow(k);
ypk=n./k(2:2:end); %Years per cycle
plot([k;k],[0*pow;pow],'c-',k,pow,'b.',…
'linewidth',2,'markersize',16)
axis([0 max(k)+1 0 pmax]);
set(gca,'xtick',k(2:2:end));
xticklabels=sprintf('%5.1f|',ypk);
set(gca,'xticklabel',xticklabels);
xlabel('years/cycle');ylabel('power')
title('Periodogram')
Sunspots
144.0 72.0 48.0 36.0 28.8 24.0 20.6 18.0 16.0 14.4 13.1 12.0 11.1 10.3 9.6 9.0 8.5 8.00
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2 x 10
7
years/cycle
po
we
r
Periodogram
原来如此
Fast Finite Fourier Transform
??
?
?? ???
1
0
11 1,.,,,0,
n
j
j
jk
k nkyY ?
每个 k需要 n个乘 n个加，总共 2n× n次浮点运算
（不包括矩阵的生成）
给出你计算机的频率，算一下需要多长时间？
当代的 FFT,1965
Jonh Tukey(Princeton University)
John Cooley(IBM Research)
)()l o g( 22 nOnnO ?
nn 2lo g3
Fastest Fourier Transform in
the West
试一下，你的机器运行 1024*1024的 fft运算要多长时间
y=rand(1024); Y=fft(y); plot(y,'.');
figure; plot(Y,'.');
Fastest Fourier Transform in the West
http://www.fftw.org MIT
Fast Finite Fourier Transform
???? ? s i nc o s/2 ie nin ???? ?
nn ?? ?22
??
???
?
?
?
?
?
?
??
?
?
??
??
???
12/
0
22
2
12/
0
12
2
odd
1
e v e n
1
1
0
11
n
j
j
jkk
n
j
j
jk
j
j
jk
j
j
jk
n
j
j
jk
k
yy
yyyY
???
???
如果 n是偶数且 k≤n/2-1,把整个求和分成两个部分
Fast Finite Fourier Transform
n=length(y);
omega=exp(-2*pi*i/n); k=(0:n/2-1)';
w=omega.^k; u=fft(y(1:2:n-1));
v=w.*fft(y(2:2:n));
fft(y)=[u+v; u-v];
如果 n是 2的 p次，这个过程可以重复 p次，每次的计算
量是 O(n),总的计算量是
)l o g()( 2 nnOnpO ?
ffttx.m
Fast Finite Fourier Transform
ffttx.m
原型：
fftmatrixgui.m
图的稀疏性，影响 FFT算法的速度
F=fft(eye(n,n)); plot(fft(eye(n,n)))
Fourier Transforms and Series
??
?
?? ???
1
0
11 1,.,,,0,
n
j
j
jk
k nkyY ?
1,...,0,1
1
0
11 ??? ?
?
?
?? njYny
n
j
jk
kj ?
? ??? ?? dtetfF ti??? 2)()(
? ???? ?? ?? deFtf ti2)()(
如果 f是周期函数,f(t+L)=f(t)
?? ? ?
?
???
?? 2/
2/
2/2 )(1,)( L
L
i j t
j
j
Li j t
j dtetfLcectf
??
