【五一创作】python序列数据处理【低通滤波】

对于海洋和气象的序列数据，可以认为其是由不同频率（波长）的信号组成的。比如30年逐小时的温度数据，有长期趋势、年际、季节和日变化，以及噪音组成。
为了研究关注的特定对象，需要将特定周期或者频率的信号提取出来，那么就需要用到滤波器。滤波器主要有低通滤波器、高通滤波器和带通滤波器3种。
低通滤波：从原始数据中，滤出低于某个频率（或某个周期）的信号，其他频率的信号剔除。
高通滤波：从原始数据中，滤出高于某个频率（或某个周期）的信号，其他频率的信号剔除。
带通滤波：从原始数据中，滤出两个频率之间的频率带，过高和过低的频率都被剔除。

这里实现一个低通滤波的例子

step1

爬取一个加拿大浮标观测的温度数据，并使用巴特沃斯滤波器（Butterworth filter），进行低通滤波。

import requests
import csv
import numpy as np
from datetime import datetime

startdate = '20110101'
enddate   = '20120101'
response = requests.get('http://lobo.satlantic.com/cgi-data/nph-data.cgi?min_date='
                           +startdate+'&max_date='+enddate+'&y=temperature')

list_data = response.text.split('\n')[1:-1] #第0个和最后一个数据不要

date = []
temp = []
for idata in list_data:
    idate, itemp = idata.split('\t')
    date.append(idate)
    temp.append(itemp)

temp = np.array(temp)
temp = temp.astype(np.float)

在这里浮标数据的网站获取为http://lobo.satlantic.com/

代码解释：这段代码使用 Python 语言，包含以下功能：

1. 导入 requests、csv 和 numpy 模块。

import requests
import csv
import numpy as np

2. 导入 datetime 模块。

from datetime import datetime

3. 定义开始日期和结束日期。

startdate = '20111118'
enddate   = '20121125'

4. 使用 requests 模块向指定 URL 发送 HTTP 请求，并获取响应数据。

response = requests.get('http://lobo.satlantic.com/cgi-data/nph-data.cgi?min_date='
                           +startdate+'&max_date='+enddate+'&y=temperature')

5. 将响应数据转换为字符串，并按行拆分成列表。

list_data = response.text.split('\n')[1:-1]

6. 使用循环遍历每行数据，并将日期和温度值分别存储到 date 和 temp 列表中。

date = []
temp = []
for idata in list_data:
    idate, itemp = idata.split('\t')
    date.append(idate)
    temp.append(itemp)

7. 将 temp 转换为 numpy 数组，并将其元素类型转换为 float。

temp = np.array(temp)
temp = temp.astype(np.float)

step2

然后进行低通滤波并绘图

import matplotlib.pyplot as plt
import scipy.signal as signal
import matplotlib.dates as mdates
# Make plots

# First, design the Buterworth filter
N  = 2    # Filter order
Wn = 0.01 # Cutoff frequency
B, A = signal.butter(N, Wn, output='ba')

# Second, apply the filter
tempf = signal.filtfilt(B,A, temp)


fig = plt.figure()
ax1 = fig.add_subplot(211)
plt.plot(date,temp, 'b-')
plt.plot(date,tempf, 'r-',linewidth=2)
plt.ylabel("Temperature (oC)")
plt.legend(['Original','Filtered'])
plt.title("Temperature from LOBO (Halifax, Canada)")

plt.xticks(year_hour_cal[0:9],[],rotation=45)
plt.xlim([0,6784])

ax1 = fig.add_subplot(212)
plt.plot(date,temp-tempf, 'b-')
plt.ylabel("Temperature (oC)")
plt.xlabel("Date")
plt.legend(['Residuals'])

plt.xticks(year_hour_cal[0:9],['Jan', 'Feb', 'Mar', 'Apr', 'May','Jun','Jul','Aug','Sep'],rotation=45)
plt.xlim([0,6784])

plt.savefig('tem_signal_filtering_plot.png', bbox_inches='tight')
plt.show()

1. 导入 scipy.signal 和 matplotlib.pyplot 模块。

import scipy.signal as signal
import matplotlib.pyplot as plt

2. 设计 Butterworth 滤波器。

N  = 2    # 滤波器阶数
Wn = 0.01 # 截止频率
B, A = signal.butter(N, Wn, output='ba')

3. 应用滤波器。

tempf = signal.filtfilt(B,A, temp)

4. 绘制图形。

fig = plt.figure()
ax1 = fig.add_subplot(211)
plt.plot(date,temp, 'b-')
plt.plot(date,tempf, 'r-',linewidth=2)
plt.ylabel("Temperature (oC)")
plt.legend(['Original','Filtered'])
plt.title("Temperature from LOBO (Halifax, Canada)")
ax1.axes.get_xaxis().set_visible(False)

ax1 = fig.add_subplot(212)
plt.plot(date,temp-tempf, 'b-')
plt.ylabel("Temperature (oC)")
plt.xlabel("Date")
plt.legend(['Residuals'])

plt.savefig(fname="./demo_1.png", dpi=300)
plt.show()

在绘图过程中，使用 fig = plt.figure() 创建一个新的图形对象，并使用 fig.add_subplot(211) 和 fig.add_subplot(212) 在图形中添加两个子图。使用 plt.plot() 绘制原始温度数据和经过滤波器处理后的温度数据，并使用 plt.legend() 添加图例。使用 plt.title() 和 plt.xlabel() 添加主标题和 x 轴标签。最后使用 plt.savefig() 将图形保存为 PNG 格式文件，并使用 plt.show() 在屏幕上显示图形。

Reference

第二篇为博主之前写的学习文章;
OceanPython+Signal filtering (Butterworth filter)
Python oceanPython项目 对温度的filtering（相当于滑动平均）【√】
python序列数据处理，低通滤波