采样率越高越好

如果你的目的是寻找频率特征，那么对信号的采样率是越高越好。

这种好并不是无脑的好，而是有数字依据的好。

采样率越高越好

采样与奈奎斯特频率

信号处理的基本理论表明，采样会导致信号的数字功率谱发生变化，采样频率的一半称为奈奎斯特频率。只有奈氏频率之下的频率成分，才能在经过正、反傅里叶变换后“复原”回来。

但这是只是“可以”，而不是“必定”，因此这是一个必要条件。

事实上，并没有人保证，经过采样和频率变换的信号，一定能“复原”回来。

本文就通过计算模拟的方式，说明一个问题，那就是

即使在奈氏频率高于感兴趣频带最高频率的情况下，采样率过低也会导致信号失真。

因此，我认为，采样频率的追求越高越好，没有上限。

实验模拟

首先，生成一个随机信号，我们假设采样频率是 1000 Hz

原信号

之后计算它的频谱，不仅如此，还对信号进行降采样到 20、50、100 Hz，它们的频谱图如下。我们有理由认为蓝色线是真实曲线，而其他颜色分别是对它进行“降采样”妥协的结果。

多降采样的频谱特征

接下来，我们通过对应频率相减的方式量化降采样带来的绝对误差，可以看到，降采样越多，则频率误差越大。表现在频率特征上，就是会多出很多莫名其妙的峰值

绝对误差

最后，我们将减值除以能量值，得到的是“比例误差”，即真实信号失真的比例有多大。可以看到，误差相当大。

比例误差

比较代码

为了可追溯起见，本文的代码如下

import sklearn.preprocessing

import numpy as np
import pandas as pd
import plotly.express as px
import plotly.subplots as subplots
import plotly.graph_objects as go

from tqdm.auto import tqdm
from noise import pnoise1

# Generate data
def generate_data(repeat=50, num=1000):

    data = []
    
    seeds = [e for e in range(num)]
    np.random.shuffle(seeds)
    seeds = seeds[:repeat]

    for seed in tqdm(seeds):
        octaves = 0 * 4 + 1 * np.random.randint(5, 7)
        persistence = 0 * 0.8 + 1 * np.random.uniform(0.7, 0.9)
        s = np.array([
            pnoise1(e, persistence=persistence, base=seed, octaves=octaves)
            for e in np.linspace(0, 1, num)
        ])
        data.append(s)

    data = np.array(data)

    data = sklearn.preprocessing.minmax_scale(data, axis=1)
    for e in data:
        e -= np.mean(e)

    print(data.shape)
    return data

data = generate_data(repeat=500)

px.imshow(data).show()

# Toolbox of comput frequency power different toolbox
def draw_diff(df):
    a = df.query('type == "Sample at 1000.0Hz"')
    b = df.query('type != "Sample at 1000.0Hz"')
    c = pd.merge(b, a, on='freq')
    c['diff'] = np.abs(c['engy_x'] - c['engy_y'])
    c['diff_ratio'] = c['diff'] / c['engy_y']
    c.index = range(len(c))
    
    c.loc[0, 'type_x']='Sample at 1000.0Hz'
    
    px.line(c, x='freq', y='diff', color='type_x', title='Diff', log_y=True).show()
    px.line(c, x='freq', y='diff_ratio', color='type_x', range_y=[0, 1], title='Diff Ratio').show()

# Compute downsample
sample_rate = 1000

def mk_downsample(downsample, signal, offset=0, sample_rate=sample_rate):
    if offset > downsample:
        offset = downsample - 1
        
    y = signal.copy()
    y = y[offset::downsample]
    fmax = sample_rate / downsample
    n = len(y)

    fft = np.fft.fftshift(np.fft.fft(y))
    freqs = np.fft.fftshift(np.fft.fftfreq(n)) * fmax

    df = pd.DataFrame()
    df['freq'] = freqs
    df['engy'] = np.abs(fft) / n
    df['type'] = 'Sample at {}Hz'.format(fmax)
    return df

# Try once
dfs = []

signal = np.mean(data, axis=0)
for down_sample in tqdm([1, 10, 20, 50], 'Simulation'):
    dfs.append(mk_downsample(down_sample, signal))

df = pd.concat(dfs)

px.line(df, x='freq', y='engy', color='type', log_y=True).show()
draw_diff(df)

# Try twice
dfs = []

signal = np.mean(data, axis=0)
for down_sample in tqdm([1, 10, 20, 50], 'Simulation'):
    dfs.append(mk_downsample(down_sample, signal, offset=5))

df = pd.concat(dfs)

px.line(df, x='freq', y='engy', color='type', log_y=True, ).show()
draw_diff(df)

# Origin signal
px.line(signal).show()