- private ArrayList m_TimeSeries = new ArrayList(); //存储时域信号的数据序列
- private ArrayList m_ZoomFFTSeries = new ArrayList(); //存储细化谱的数据序列
- //tempTr、tempTi暂存kbfft()的入口数据;tempFr、tempFi存储kbfft()返回的实部和虚部
- private ArrayList tempTr = new ArrayList(), tempTi = new ArrayList(), tempFr = new ArrayList(), tempFi = new ArrayList();
- private double m_Xstep, m_ZoomXmax, m_ZoomXmin;
- public ArrayList TimeSeries
- {
- get{ return m_TimeSeries; }
- set{ m_TimeSeries = value; }
- }
- public ArrayList ZoomFFTSeries
- {
- get { return m_ZoomFFTSeries; }
- set { m_ZoomFFTSeries = value; }
- }
- public double Xstep
- {
- get { return m_Xstep; }
- set { m_Xstep = value; }
- }
- public double ZoomXmax
- {
- get { return m_ZoomXmax; }
- set { m_ZoomXmax = value; }
- }
- public double ZoomXmin
- {
- get { return m_ZoomXmin; }
- set { m_ZoomXmin = value; }
- }
-
- #region 细化谱序列
- /// <summary>
- /// 通过时域信号的数据序列,计算细化谱数据序列。fin为分析中心频率,iN为做谱点数,scale细化倍数,iL平均段数,iM为滤波器半阶数
- /// </summary>
- public void GetZoomFFTSeries(double fin, int scale, int iN, int iL, int iM)//
- {
- m_ZoomFFTSeries.Clear();
- double[] w = new double[iM];
- for (int i = 0; i < iM; i++)
- {
- w[i] = 0.5 + 0.5 * Math.Cos(Math.PI * (i+1) / iM);//Hanning窗
- }
- double fl, fh;//最小截止频率、最大截止频率
- //fl = (fin - m_SampleHz / (4.0 * scale)) > (-m_SampleHz / 2.2) ? (fin - m_SampleHz / (4 * scale)) : -m_SampleHz / 2.2;
- //fh = (fin + m_SampleHz / (4.0 * scale)) < (m_SampleHz / 2.2) ? (fin + m_SampleHz / (4 * scale)) : m_SampleHz / 2.2;
- fl = fin - m_SampleHz / (4.0 * scale); fh = fin + m_SampleHz / (4.0 * scale);
- double yf = scale * fl;//移频量
- double df = m_SampleHz / scale / iN;//
- double[] xz = new double[iN / 2];
- double wl = 2 * Math.PI * fl / m_SampleHz, wh = 2 * Math.PI * fh / m_SampleHz;
- double[] hr = new double[iM + 1], hi = new double[iM + 1];
- hr[0] = (wl - wh) / Math.PI; hi[0] = 0;
- for (int i = 1; i <= iM; i++)
- {
- hr[i] = (Math.Sin(wl * i) - Math.Sin(wh * i)) / (Math.PI * i) * w[i - 1];
- hi[i] = (Math.Cos(wl * i) - Math.Cos(wh * i)) / (Math.PI * i) * w[i - 1];
- }
- double[] w1 = new double[iN];
- for (int i = 0; i < iN; i++)
- {
- w1[i] = 0.5 - 0.5 * Math.Cos(2 * Math.PI * i / iN);
- }
- double[] xrz = new double[iN], xiz = new double[iN],xzt=new double[2*iN];
- double[] x = new double[m_TimeSeries.Count];
- for (int i = 0; i < m_TimeSeries.Count; i++)
- {
- x[i] = Convert.ToDouble(m_TimeSeries[i]);
- }
- double factor = -2 * Math.PI * yf / m_SampleHz;
- for (int i = 1; i <= iL; i++)
- {
- for (int k = 1; k <= iN; k++)
- {
- int kk = (k - 1) * scale + iM + (i - 1) * iN;
- double sumR = 0, sumI = 0;
- for (int j = 1; j <= iM; j++)
- {
- sumR += hr[j] * (x[kk + j] + x[kk - j]);
- sumI += hi[j] * (x[kk + j] - x[kk - j]);
- }
- xrz[k-1]=(x[kk] * hr[0] + sumR);
- xiz[k-1]=(x[kk] * hi[0] + sumI);
- }
- double sumXZTr = 0, sumXZTi = 0;
- for (int k = 0; k < 2*iN; k+=2)
- {
- xzt[k] = xrz[k / 2] * Math.Cos(k * factor / 2) - xiz[k / 2] * Math.Sin(k * factor / 2);//偶数位为实部,奇数位为虚部
- xzt[k + 1] = xiz[k / 2] * Math.Cos(k * factor / 2) + xrz[k / 2] * Math.Sin(k * factor / 2);
- xzt[k] = xzt[k] * w1[k / 2];
- xzt[k + 1] = xzt[k + 1] * w1[k / 2];
- sumXZTr += xzt[k]; sumXZTi += xzt[k + 1];
- }
- tempTr.Clear(); tempTi.Clear();
- for (int k = 0; k < 2 * iN; k += 2)
- {
- tempTr.Add(xzt[k] - sumXZTr / iN);
- tempTi.Add(xzt[k + 1] - sumXZTi / iN);
- }
- kbfft(tempTr.Count, 0, 0);
- for (int k = 0; k < iN / 2; k++)
- {
- double temR = Convert.ToDouble(tempFr[k]), temI = Convert.ToDouble(tempFi[k]);
- xz[k] += Math.Pow(Math.Sqrt(temI * temI + temR * temR) / iN * 2, 2);
- }
- }
- for (int i = 0; i < iN / 2; i++)
- {
- m_ZoomFFTSeries.Add(Math.Sqrt(xz[i] / iL));
- }
- m_ZoomXmin = fl; m_ZoomXmax = fh;
- m_Xstep = (fh - fl) / iN * 2;
- }
-
- /// <summary>
- /// GetZoomFFTSeries(double fin, int scale, int iN, int iL, int iM)的重构方法。
- /// </summary>
- public void GetZoomFFTSeries()
- {
- this.GetZoomFFTSeries(500, 5, 1024, 1, m_TimeSeries.Count/16);
- }
- #endregion
- #region 快速傅里叶变换方法
- /***************************************************************************
- * 入口参数:
- * l: l = 0, 傅立叶变换; l = 1, 逆傅立叶变换
- * il: il = 0,不计算傅立叶变换或逆变换模和幅角;il = 1,计算模和幅角
- * n: 输入的点数,为偶数,一般为32,64,128,...,1024等
- * k: 满足n=2^k(k>0),实质上k是n个采样数据可以分解为偶次幂和奇次幂的次数
- * pr[]: l=0时,存放N点采样数据的实部
- * l=1时, 存放傅立叶变换的N个实部
- * pi[]: l=0时,存放N点采样数据的虚部
- * l=1时, 存放傅立叶变换的N个虚部
- *
- * 出口参数:
- * fr[]: l=0, 返回傅立叶变换的实部
- * l=1, 返回逆傅立叶变换的实部
- * fi[]: l=0, 返回傅立叶变换的虚部
- * l=1, 返回逆傅立叶变换的虚部
- * pr[]: il = 1,l = 0 时,返回傅立叶变换的模
- * il = 1,l = 1 时,返回逆傅立叶变换的模
- * pi[]: il = 1,l = 0 时,返回傅立叶变换的辐角
- * il = 1,l = 1 时,返回逆傅立叶变换的辐角
- ***************************************************************************/
- /// <summary>
- /// 快速傅里叶变换方法。
- /// </summary>
- void kbfft(int n,int l, int il)
- {
- int k = (int)Math.Log(n, 2);
- double[] pr = new double[n], pi = new double[n], fr = new double[n], fi = new double[n];
- for (int i = 0; i < n; i++)
- {
- if (i < tempTr.Count)
- pr[i] = Convert.ToDouble(tempTr[i]);
- else
- pr[i] = 0;
- if (i < tempTi.Count)
- pi[i] = Convert.ToDouble(tempTi[i]);
- else
- pi[i] = 0;
- }
- int m, is1, j, nv;
- double p, q, s, vr, vi, poddr, poddi;
- //排序
- for (int it = 0; it <= n - 1; it++)
- {
- m = it; is1 = 0;
- for (int i = 0; i <= k - 1; i++)
- {
- j = m / 2; is1 = 2 * is1 + (m - 2 * j); m = j;
- fr[it] = pr[is1]; fi[it] = pi[is1];
- }
- }
-
- //蝶形运算
- pr[0] = 1.0; pi[0] = 0.0;
- p = Math.PI * 2 / (1.0 * n);
- pr[1] = Math.Cos(p); pi[1] = -Math.Sin(p);
- if (l != 0) pi[1] = -pi[1];
- for (int i = 2; i <= n - 1; i++)
- {
- p = pr[i - 1] * pr[1]; q = pi[i - 1] * pi[1];
- s = (pr[i - 1] + pi[i - 1]) * (pr[1] + pi[1]);
- pr[i] = p - q; pi[i] = s - p - q;
- }
- for (int it = 0; it <= n - 2; it = it + 2)
- {
- vr = fr[it]; vi = fi[it];
- fr[it] = vr + fr[it + 1]; fi[it] = vi + fi[it + 1];
- fr[it + 1] = vr - fr[it + 1]; fi[it + 1] = vi - fi[it + 1];
- }
- m = n / 2; nv = 2;
- for (int l0 = k - 2; l0 >= 0; l0--)
- {
- m = m / 2; nv = 2 * nv;
- for (int it = 0; it <= (m - 1) * nv; it = it + nv)
- for ( j = 0; j <= (nv / 2) - 1; j++)
- {
- p = pr[m * j] * fr[it + j + nv / 2];
- q = pi[m * j] * fi[it + j + nv / 2];
- s = pr[m * j] + pi[m * j];
- s = s * (fr[it + j + nv / 2] + fi[it + j + nv / 2]);
- poddr = p - q; poddi = s - p - q;
- fr[it + j + nv / 2] = fr[it + j] - poddr;
- fi[it + j + nv / 2] = fi[it + j] - poddi;
- fr[it + j] = fr[it + j] + poddr;
- fi[it + j] = fi[it + j] + poddi;
- }
- }
- if (l != 0)
- for (int i = 0; i <= n - 1; i++)
- {
- fr[i] = fr[i] / (1.0 * n);
- fi[i] = fi[i] / (1.0 * n);
- }
- if (il != 0)
- for (int i = 0; i <= n - 1; i++)
- {
- pr[i] = Math.Sqrt(fr[i] * fr[i] + fi[i] * fi[i]);
- pr[i] = (pr[i] / (n / 2)); //各次谐波幅值,其中pr[1]为基波幅值
- if (Math.Abs(fr[i]) < 0.000001 * Math.Abs(fi[i]))
- {
- if ((fi[i] * fr[i]) > 0) pi[i] = 90.0;
- else pi[i] = -90.0;
- }
- else
- pi[i] = Math.Atan(fi[i] / fr[i]) * 360.0 / (2*Math.PI);
- }
- tempFi.Clear(); tempFr.Clear();
- for (int i = 0; i < n; i++)
- {
- tempFr.Add(fr[i]);
- tempFi.Add(fi[i]);
- }
- }
- #endregion
复制代码 这是我写的一个类(其中一部分)。。。希望对大家有所帮助 |