声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

楼主: yangzj

[FFT] [原创]zoomfft

[复制链接]
发表于 2007-4-12 22:43 | 显示全部楼层

明白了,谢谢
回复 支持 反对
分享到:

使用道具 举报

发表于 2007-6-18 11:09 | 显示全部楼层
现在我是一点也看不懂,不过我相信以后就会懂了。
发表于 2008-2-27 10:33 | 显示全部楼层
hr(1)=(wl-wh)/pi;
hr(2:M+1)=(sin(wl*k)-sin(wh*k))./(pi*k).*w;
hi(1)=0;
hi(2:M+1)=(cos(wl*k)-cos(wh*k))./(pi*k).*w;

程序是不是有误,对照丁康的论文,缺少了k=1时的值,却多了k=M+1时的值?
发表于 2008-2-27 10:56 | 显示全部楼层

晕了,看明白了,不好意思

晕了,看明白了,不好意思
发表于 2009-6-17 03:58 | 显示全部楼层

zoomfft函数不能调用

怎么会出现“??? Undefined command/function 'zoomfft'.“,高手教我!
发表于 2009-6-17 08:17 | 显示全部楼层

回复 50楼 JLBhaidao 的帖子

使用which zoomfft检查下
不是无此函数就是路径没设好
发表于 2010-2-10 20:10 | 显示全部楼层
弱弱的问一下
for k=1:N
        kk=(k-1)*D+M+(i-1)*N;
        xrz(k)=x(kk+1)*hr(1)+sum(hr(2:M+1).*(x(kk+2:kk+M+1)+x(kk:-1:kk-M+1)));
        xiz(k)=x(kk+1)*hi(1)+sum(hi(2:M+1).*(x(kk+2:kk+M+1)-x(kk:-1:kk-M+1)));
    end
这段程序是在进行滤波和抽取吗?
发表于 2010-2-11 14:17 | 显示全部楼层
学习了:lol
发表于 2010-4-9 22:30 | 显示全部楼层

回复 楼主 yangzj 的帖子

请教给位仁兄,我运行这个程序时,会出现以下错误信号,不知为何?请指教,多谢!
??? function [f xz]=ZoomFFT(x,fs,N,fe,D,L,M)
    |
Error: Function definitions are not permitted at the prompt or in scripts.
发表于 2010-8-4 21:05 | 显示全部楼层
学习了:@L 还是不太懂
发表于 2010-8-24 00:07 | 显示全部楼层
谢谢楼主,拷下来学习了
发表于 2010-8-26 11:30 | 显示全部楼层
回复 shanghai 的帖子


   
发表于 2010-11-30 16:44 | 显示全部楼层
发表于 2011-4-16 13:18 | 显示全部楼层
纠结这个细化谱好几天了,也找了网上和书上其它的matlab源程序,但改写C#始终不对。。。今天用楼主的这个就对了~~非常感谢。。
发表于 2011-4-16 13:23 | 显示全部楼层
  1.     private ArrayList m_TimeSeries = new ArrayList();           //存储时域信号的数据序列
  2.     private ArrayList m_ZoomFFTSeries = new ArrayList();        //存储细化谱的数据序列
  3.     //tempTr、tempTi暂存kbfft()的入口数据;tempFr、tempFi存储kbfft()返回的实部和虚部
  4.     private ArrayList tempTr = new ArrayList(), tempTi = new ArrayList(), tempFr = new ArrayList(), tempFi = new ArrayList();
  5.     private double m_Xstep, m_ZoomXmax, m_ZoomXmin;

  6.     public ArrayList TimeSeries
  7.     {
  8.         get{ return m_TimeSeries; }
  9.         set{ m_TimeSeries = value; }
  10.     }
  11.     public ArrayList ZoomFFTSeries
  12.     {
  13.         get { return m_ZoomFFTSeries; }
  14.         set { m_ZoomFFTSeries = value; }
  15.     }
  16.     public double Xstep
  17.     {
  18.         get { return m_Xstep; }
  19.         set { m_Xstep = value; }
  20.     }
  21.     public double ZoomXmax
  22.     {
  23.         get { return m_ZoomXmax; }
  24.         set { m_ZoomXmax = value; }
  25.     }
  26.     public double ZoomXmin
  27.     {
  28.         get { return m_ZoomXmin; }
  29.         set { m_ZoomXmin = value; }
  30.     }


  31.     #region 细化谱序列
  32.     /// <summary>
  33.     /// 通过时域信号的数据序列,计算细化谱数据序列。fin为分析中心频率,iN为做谱点数,scale细化倍数,iL平均段数,iM为滤波器半阶数
  34.     /// </summary>
  35.     public void GetZoomFFTSeries(double fin, int scale, int iN, int iL, int iM)//
  36.     {
  37.         m_ZoomFFTSeries.Clear();
  38.         double[] w = new double[iM];
  39.         for (int i = 0; i < iM; i++)
  40.         {
  41.             w[i] = 0.5 + 0.5 * Math.Cos(Math.PI * (i+1) / iM);//Hanning窗
  42.         }
  43.         double fl, fh;//最小截止频率、最大截止频率
  44.         //fl = (fin - m_SampleHz / (4.0 * scale)) > (-m_SampleHz / 2.2) ? (fin - m_SampleHz / (4 * scale)) : -m_SampleHz / 2.2;
  45.         //fh = (fin + m_SampleHz / (4.0 * scale)) < (m_SampleHz / 2.2) ? (fin + m_SampleHz / (4 * scale)) : m_SampleHz / 2.2;
  46.         fl = fin - m_SampleHz / (4.0 * scale); fh = fin + m_SampleHz / (4.0 * scale);

  47.         double yf = scale * fl;//移频量
  48.         double df = m_SampleHz / scale / iN;//

  49.         double[] xz = new double[iN / 2];
  50.         double wl = 2 * Math.PI * fl / m_SampleHz, wh = 2 * Math.PI * fh / m_SampleHz;
  51.         double[] hr = new double[iM + 1], hi = new double[iM + 1];
  52.         hr[0] = (wl - wh) / Math.PI; hi[0] = 0;
  53.         for (int i = 1; i <= iM; i++)
  54.         {
  55.             hr[i] = (Math.Sin(wl * i) - Math.Sin(wh * i)) / (Math.PI * i) * w[i - 1];
  56.             hi[i] = (Math.Cos(wl * i) - Math.Cos(wh * i)) / (Math.PI * i) * w[i - 1];
  57.         }
  58.         double[] w1 = new double[iN];
  59.         for (int i = 0; i < iN; i++)
  60.         {
  61.             w1[i] = 0.5 - 0.5 * Math.Cos(2 * Math.PI * i / iN);
  62.         }

  63.         double[] xrz = new double[iN], xiz = new double[iN],xzt=new double[2*iN];
  64.         double[] x = new double[m_TimeSeries.Count];
  65.         for (int i = 0; i < m_TimeSeries.Count; i++)
  66.         {
  67.             x[i] = Convert.ToDouble(m_TimeSeries[i]);
  68.         }
  69.         double factor = -2 * Math.PI * yf / m_SampleHz;
  70.         for (int i = 1; i <= iL; i++)
  71.         {
  72.             for (int k = 1; k <= iN; k++)
  73.             {
  74.                 int kk = (k - 1) * scale + iM + (i - 1) * iN;
  75.                 double sumR = 0, sumI = 0;
  76.                 for (int j = 1; j <= iM; j++)
  77.                 {
  78.                     sumR += hr[j] * (x[kk + j] + x[kk - j]);
  79.                     sumI += hi[j] * (x[kk + j] - x[kk - j]);
  80.                 }
  81.                 xrz[k-1]=(x[kk] * hr[0] + sumR);
  82.                 xiz[k-1]=(x[kk] * hi[0] + sumI);
  83.             }
  84.             double sumXZTr = 0, sumXZTi = 0;
  85.             for (int k = 0; k < 2*iN; k+=2)
  86.             {
  87.                 xzt[k] = xrz[k / 2] * Math.Cos(k * factor / 2) - xiz[k / 2] * Math.Sin(k * factor / 2);//偶数位为实部,奇数位为虚部
  88.                 xzt[k + 1] = xiz[k / 2] * Math.Cos(k * factor / 2) + xrz[k / 2] * Math.Sin(k * factor / 2);
  89.                 xzt[k] = xzt[k] * w1[k / 2];
  90.                 xzt[k + 1] = xzt[k + 1] * w1[k / 2];
  91.                 sumXZTr += xzt[k]; sumXZTi += xzt[k + 1];
  92.             }
  93.             tempTr.Clear(); tempTi.Clear();
  94.             for (int k = 0; k < 2 * iN; k += 2)
  95.             {
  96.                 tempTr.Add(xzt[k] - sumXZTr / iN);
  97.                 tempTi.Add(xzt[k + 1] - sumXZTi / iN);
  98.             }
  99.             kbfft(tempTr.Count, 0, 0);
  100.             for (int k = 0; k < iN / 2; k++)
  101.             {
  102.                 double temR = Convert.ToDouble(tempFr[k]), temI = Convert.ToDouble(tempFi[k]);
  103.                 xz[k] += Math.Pow(Math.Sqrt(temI * temI + temR * temR) / iN * 2, 2);
  104.             }
  105.         }
  106.         for (int i = 0; i < iN / 2; i++)
  107.         {
  108.             m_ZoomFFTSeries.Add(Math.Sqrt(xz[i] / iL));
  109.         }
  110.         m_ZoomXmin = fl; m_ZoomXmax = fh;
  111.         m_Xstep = (fh - fl) / iN * 2;
  112.     }
  113.    
  114.     /// <summary>
  115.     /// GetZoomFFTSeries(double fin, int scale, int iN, int iL, int iM)的重构方法。
  116.     /// </summary>
  117.     public void GetZoomFFTSeries()
  118.     {
  119.         this.GetZoomFFTSeries(500, 5, 1024, 1, m_TimeSeries.Count/16);
  120.     }
  121.     #endregion

  122.     #region 快速傅里叶变换方法
  123.     /***************************************************************************
  124.      * 入口参数:
  125.      *  l: l = 0, 傅立叶变换; l = 1, 逆傅立叶变换
  126.      * il: il = 0,不计算傅立叶变换或逆变换模和幅角;il = 1,计算模和幅角
  127.      * n: 输入的点数,为偶数,一般为32,64,128,...,1024等
  128.      * k: 满足n=2^k(k>0),实质上k是n个采样数据可以分解为偶次幂和奇次幂的次数
  129.      * pr[]: l=0时,存放N点采样数据的实部
  130.      * l=1时, 存放傅立叶变换的N个实部
  131.      * pi[]: l=0时,存放N点采样数据的虚部
  132.      * l=1时, 存放傅立叶变换的N个虚部
  133.      *
  134.      * 出口参数:
  135.      * fr[]: l=0, 返回傅立叶变换的实部
  136.      * l=1, 返回逆傅立叶变换的实部
  137.      * fi[]: l=0, 返回傅立叶变换的虚部
  138.      * l=1, 返回逆傅立叶变换的虚部
  139.      * pr[]: il = 1,l = 0 时,返回傅立叶变换的模
  140.      * il = 1,l = 1 时,返回逆傅立叶变换的模
  141.      * pi[]: il = 1,l = 0 时,返回傅立叶变换的辐角
  142.      * il = 1,l = 1 时,返回逆傅立叶变换的辐角
  143.      ***************************************************************************/
  144.     /// <summary>
  145.     /// 快速傅里叶变换方法。
  146.     /// </summary>
  147.     void kbfft(int n,int l, int il)
  148.     {
  149.         int k = (int)Math.Log(n, 2);
  150.         double[] pr = new double[n], pi = new double[n], fr = new double[n], fi = new double[n];
  151.         for (int i = 0; i < n; i++)
  152.         {
  153.             if (i < tempTr.Count)
  154.                 pr[i] = Convert.ToDouble(tempTr[i]);
  155.             else
  156.                 pr[i] = 0;
  157.             if (i < tempTi.Count)
  158.                 pi[i] = Convert.ToDouble(tempTi[i]);
  159.             else
  160.                 pi[i] = 0;
  161.         }

  162.         int  m, is1, j, nv;
  163.         double p, q, s, vr, vi, poddr, poddi;

  164.         //排序
  165.         for (int it = 0; it <= n - 1; it++)
  166.         {
  167.             m = it; is1 = 0;
  168.             for (int i = 0; i <= k - 1; i++)
  169.             {
  170.                 j = m / 2; is1 = 2 * is1 + (m - 2 * j); m = j;
  171.                 fr[it] = pr[is1]; fi[it] = pi[is1];
  172.             }
  173.         }
  174.         
  175.         //蝶形运算
  176.         pr[0] = 1.0; pi[0] = 0.0;
  177.         p = Math.PI * 2 / (1.0 * n);
  178.         pr[1] = Math.Cos(p); pi[1] = -Math.Sin(p);
  179.         if (l != 0) pi[1] = -pi[1];
  180.         for (int i = 2; i <= n - 1; i++)
  181.         {
  182.             p = pr[i - 1] * pr[1]; q = pi[i - 1] * pi[1];
  183.             s = (pr[i - 1] + pi[i - 1]) * (pr[1] + pi[1]);
  184.             pr[i] = p - q; pi[i] = s - p - q;
  185.         }
  186.         for (int it = 0; it <= n - 2; it = it + 2)
  187.         {
  188.             vr = fr[it]; vi = fi[it];
  189.             fr[it] = vr + fr[it + 1]; fi[it] = vi + fi[it + 1];
  190.             fr[it + 1] = vr - fr[it + 1]; fi[it + 1] = vi - fi[it + 1];
  191.         }
  192.         m = n / 2; nv = 2;
  193.         for (int l0 = k - 2; l0 >= 0; l0--)
  194.         {
  195.             m = m / 2; nv = 2 * nv;
  196.             for (int it = 0; it <= (m - 1) * nv; it = it + nv)
  197.                 for ( j = 0; j <= (nv / 2) - 1; j++)
  198.                 {
  199.                     p = pr[m * j] * fr[it + j + nv / 2];
  200.                     q = pi[m * j] * fi[it + j + nv / 2];
  201.                     s = pr[m * j] + pi[m * j];
  202.                     s = s * (fr[it + j + nv / 2] + fi[it + j + nv / 2]);
  203.                     poddr = p - q; poddi = s - p - q;
  204.                     fr[it + j + nv / 2] = fr[it + j] - poddr;
  205.                     fi[it + j + nv / 2] = fi[it + j] - poddi;
  206.                     fr[it + j] = fr[it + j] + poddr;
  207.                     fi[it + j] = fi[it + j] + poddi;
  208.                 }
  209.         }
  210.         if (l != 0)
  211.             for (int i = 0; i <= n - 1; i++)
  212.             {
  213.                 fr[i] = fr[i] / (1.0 * n);
  214.                 fi[i] = fi[i] / (1.0 * n);
  215.             }
  216.         if (il != 0)
  217.             for (int i = 0; i <= n - 1; i++)
  218.             {
  219.                 pr[i] = Math.Sqrt(fr[i] * fr[i] + fi[i] * fi[i]);
  220.                 pr[i] = (pr[i] / (n / 2)); //各次谐波幅值,其中pr[1]为基波幅值
  221.                 if (Math.Abs(fr[i]) < 0.000001 * Math.Abs(fi[i]))
  222.                 {
  223.                     if ((fi[i] * fr[i]) > 0) pi[i] = 90.0;
  224.                     else pi[i] = -90.0;
  225.                 }
  226.                 else
  227.                     pi[i] = Math.Atan(fi[i] / fr[i]) * 360.0 / (2*Math.PI);
  228.             }


  229.         tempFi.Clear(); tempFr.Clear();
  230.         for (int i = 0; i < n; i++)
  231.         {
  232.             tempFr.Add(fr[i]);
  233.             tempFi.Add(fi[i]);
  234.         }
  235.     }
  236.     #endregion
复制代码
这是我写的一个类(其中一部分)。。。希望对大家有所帮助
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-6 07:45 , Processed in 0.071993 second(s), 16 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表