谐波分析

1.txt中只有一列数据，共240个

#include "iostream"
#include "fstream"
#include "math.h"
#include "algorithm"
#include "vector"
using namespace std;

#define pi 3.141592653

struct V
{
    double a, b, I;
    int k;

    V(double aa, double bb, double ii, int kk)
    {
        a = aa;
        b = bb;
        I = ii;
        k = kk;
    }

    bool operator < (struct V &one)
    {
        return one.I < I;
    }
};

int main()
{
    ifstream fin;
    fin.open("1.txt");
    int n = 0;
    vector<double> x;
    double a0 = 0;
    int c;
    while(!fin.eof())
    {
        fin >> c;
        a0 += c;
        x.push_back((double)c);
        n++;
    }
    a0 /= n;
    int k = n / 2;   //波的个数
    vector<V> r;
    int i, j;
    double Isum = 0;
    for(j=1; j<=k; j++)
    {
        double a, b, I;
        a = b = I = 0;
        for(int t=0; t<n; t++)
        {
            a += x[t] * cos(2.0 * pi * j * t / (double)n);
            b += x[t] * sin(2.0 * pi * j * t / (double)n);
        }
        a = a * 2.0 / (double)n;
        b = a * 2.0 / (double)n;
        I = (a * a + b * b) / 2.0;
        V v(a, b, I, j);
        Isum += I;
        r.push_back(v);
    }

    sort(r.begin(), r.end());
    vector<V> result;
    for(i=0; i<k; i++)
    {
        r[i].I /= Isum;
        if(r[i].I > 0.071)   //如果显著
            result.push_back(r[i]);
        else
            break;
    }

    for(i=0; i<result.size(); i++)
    {
        cout << result[i].a << " " << result[i].b << " ";
        cout << result[i].k << " " << result[i].I << endl;
    }
    cout << a0 << endl;
    cout << n << endl;
    fin.close();
    return 0;
}

n = 240     %数据个数
t = 1: 1: n;
y = 275.354  - 216.137 * cos(20 * pi * t / n) - 1.801 * sin(20 * pi * t / n) ... 
             - 85.665 * cos(100 * pi * t / n) - 0.714 * sin(100 * pi * t / n);  %波方程
x = importdata('5.txt');  %导入实际数据
x = x.';    %转置
subplot(2, 1, 1);
plot(t, x);
subplot(2, 1, 2);
plot(t, y);

上面为原始数据，下面为提取出的波

参考文献
李平，卢文喜，杨忠平. 频谱分析法在吉林西部地下水动态预报中的应用[J]. 水文地质工程，2005，(4)：70-73.