上次以 Python 搭配 matplotlib 改寫張智星老師傅立葉轉換教學例子。後來逛到 Anders Andreasen 的專文，裡面有個分析太陽黑子活動週期的例子，相同的例子竟然也出現在 Mathworks 展示 Matlab FFT 用法的網頁上。既然大家那麼愛用太陽黑子，我也來攪和攪和，再次以 Python 搭配 matplotlib 改寫：
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 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 |
#!/usr/bin/python """ This demonstration uses the FFT function to analyze the variations in sunspot activity over the last 300 years. Sunspot activity is cyclical, reaching a maximum about every 11 years. Let's confirm that. Here is a plot of a quantity called the Wolfer number, which measures both number and size of sunspots. Astronomers have tabulated this number for almost 300 years. Note: This demonstration is original from both Anders Andreasen's and Mathworks' page. I translated it into Python with matplotlib. See also: - "Python for scientific use, Part II: Data analysis" by Anders Andreasen <http://linuxgazette.net/115/andreasen.html> - Using FFT in Matlab <http://www.mathworks.com/products/demos/shipping/matlab/sunspots.html> """ __author__ = "Jiang Yu-Kuan, yukuan.jiang(at)gmail.com" __date__ = "December 2006" __revision__ = "1.2" import scipy.io.array_import from pylab import * figure(figsize=(13,4.5)) subplot(121) sunspot = scipy.io.array_import.read_array('sunspots.dat') year = sunspot[:,0] wolfer = sunspot[:,1] plot(year, wolfer, "r+-") xlabel('Year') ylabel('Wolfer number') title('Sunspot data') subplot(122) Y = fft(wolfer) plot(Y.real, Y.imag, "ro") xlabel('Real Axis') ylabel('Imaginary Axis') title('Fourier Coefficients in the Complex Plane') xlim(-4000, 2000) figure(figsize=(13,4.5)) subplot(121) N = len(Y) power = abs(Y[:(N/2)])**2 Fs = 1. nyquist = Fs/2 freq = linspace(0,1,N/2)*nyquist plot(freq[1:], power[1:]) xlabel('Frequency (Cycles/Year)') ylabel('Power') title("Spectrum") xlim(0, 0.20) subplot(122) period = 1./freq plot(period[1:], power[1:]) index = find(power==max(power[1:])) plot(period[index], power[index], 'ro') text(float(period[index])+1, float(power[index])*.95, 'Period='+`float(period[index])`) xlabel('Period (Years/Cycle)') ylabel('Power') title("Periodogram") xlim(0, 40) show() |
程式會把三百年份的太陽黑子資料檔調進來分析。結果截圖如下：
0 comments:
Post a Comment