[原创]zoomfft

yixiaofan

噢
明白了，谢谢

zhaoyixu

现在我是一点也看不懂，不过我相信以后就会懂了。

xiaodinger

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时的值？

xiaodinger

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

JLBhaidao

怎么会出现“??? Undefined command/function 'zoomfft'.“，高手教我！

ChaChing

使用which zoomfft检查下
不是无此函数就是路径没设好

diandidiandi

弱弱的问一下
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
这段程序是在进行滤波和抽取吗？

无限剑制

学习了:lol

zmp007

请教给位仁兄，我运行这个程序时，会出现以下错误信号，不知为何？请指教，多谢！
??? function [f xz]=ZoomFFT(x,fs,N,fe,D,L,M)
    |
Error: Function definitions are not permitted at the prompt or in scripts.

zxk5511

学习了:@L 还是不太懂

secondye

谢谢楼主，拷下来学习了

土木年华

回复 shanghai 的帖子
牛

   

200810100659

ylhlly

纠结这个细化谱好几天了，也找了网上和书上其它的matlab源程序，但改写C#始终不对。。。今天用楼主的这个就对了~~非常感谢。。

ylhlly

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. 
   
7.     public ArrayList TimeSeries
   
8.     {
   
9.         get{ return m_TimeSeries; }
   
10.         set{ m_TimeSeries = value; }
    
11.     }
    
12.     public ArrayList ZoomFFTSeries
    
13.     {
    
14.         get { return m_ZoomFFTSeries; }
    
15.         set { m_ZoomFFTSeries = value; }
    
16.     }
    
17.     public double Xstep
    
18.     {
    
19.         get { return m_Xstep; }
    
20.         set { m_Xstep = value; }
    
21.     }
    
22.     public double ZoomXmax
    
23.     {
    
24.         get { return m_ZoomXmax; }
    
25.         set { m_ZoomXmax = value; }
    
26.     }
    
27.     public double ZoomXmin
    
28.     {
    
29.         get { return m_ZoomXmin; }
    
30.         set { m_ZoomXmin = value; }
    
31.     }
    
32. 
    
33. 
    
34.     #region 细化谱序列
    
35.     /// <summary>
    
36.     /// 通过时域信号的数据序列，计算细化谱数据序列。fin为分析中心频率，iN为做谱点数，scale细化倍数，iL平均段数，iM为滤波器半阶数
    
37.     /// </summary>
    
38.     public void GetZoomFFTSeries(double fin, int scale, int iN, int iL, int iM)//
    
39.     {
    
40.         m_ZoomFFTSeries.Clear();
    
41.         double[] w = new double[iM];
    
42.         for (int i = 0; i < iM; i++)
    
43.         {
    
44.             w[i] = 0.5 + 0.5 * Math.Cos(Math.PI * (i+1) / iM);//Hanning窗
    
45.         }
    
46.         double fl, fh;//最小截止频率、最大截止频率
    
47.         //fl = (fin - m_SampleHz / (4.0 * scale)) > (-m_SampleHz / 2.2) ? (fin - m_SampleHz / (4 * scale)) : -m_SampleHz / 2.2;
    
48.         //fh = (fin + m_SampleHz / (4.0 * scale)) < (m_SampleHz / 2.2) ? (fin + m_SampleHz / (4 * scale)) : m_SampleHz / 2.2;
    
49.         fl = fin - m_SampleHz / (4.0 * scale); fh = fin + m_SampleHz / (4.0 * scale);
    
50. 
    
51.         double yf = scale * fl;//移频量
    
52.         double df = m_SampleHz / scale / iN;//
    
53. 
    
54.         double[] xz = new double[iN / 2];
    
55.         double wl = 2 * Math.PI * fl / m_SampleHz, wh = 2 * Math.PI * fh / m_SampleHz;
    
56.         double[] hr = new double[iM + 1], hi = new double[iM + 1];
    
57.         hr[0] = (wl - wh) / Math.PI; hi[0] = 0;
    
58.         for (int i = 1; i <= iM; i++)
    
59.         {
    
60.             hr[i] = (Math.Sin(wl * i) - Math.Sin(wh * i)) / (Math.PI * i) * w[i - 1];
    
61.             hi[i] = (Math.Cos(wl * i) - Math.Cos(wh * i)) / (Math.PI * i) * w[i - 1];
    
62.         }
    
63.         double[] w1 = new double[iN];
    
64.         for (int i = 0; i < iN; i++)
    
65.         {
    
66.             w1[i] = 0.5 - 0.5 * Math.Cos(2 * Math.PI * i / iN);
    
67.         }
    
68. 
    
69.         double[] xrz = new double[iN], xiz = new double[iN],xzt=new double[2*iN];
    
70.         double[] x = new double[m_TimeSeries.Count];
    
71.         for (int i = 0; i < m_TimeSeries.Count; i++)
    
72.         {
    
73.             x[i] = Convert.ToDouble(m_TimeSeries[i]);
    
74.         }
    
75.         double factor = -2 * Math.PI * yf / m_SampleHz;
    
76.         for (int i = 1; i <= iL; i++)
    
77.         {
    
78.             for (int k = 1; k <= iN; k++)
    
79.             {
    
80.                 int kk = (k - 1) * scale + iM + (i - 1) * iN;
    
81.                 double sumR = 0, sumI = 0;
    
82.                 for (int j = 1; j <= iM; j++)
    
83.                 {
    
84.                     sumR += hr[j] * (x[kk + j] + x[kk - j]);
    
85.                     sumI += hi[j] * (x[kk + j] - x[kk - j]);
    
86.                 }
    
87.                 xrz[k-1]=(x[kk] * hr[0] + sumR);
    
88.                 xiz[k-1]=(x[kk] * hi[0] + sumI);
    
89.             }
    
90.             double sumXZTr = 0, sumXZTi = 0;
    
91.             for (int k = 0; k < 2*iN; k+=2)
    
92.             {
    
93.                 xzt[k] = xrz[k / 2] * Math.Cos(k * factor / 2) - xiz[k / 2] * Math.Sin(k * factor / 2);//偶数位为实部，奇数位为虚部
    
94.                 xzt[k + 1] = xiz[k / 2] * Math.Cos(k * factor / 2) + xrz[k / 2] * Math.Sin(k * factor / 2);
    
95.                 xzt[k] = xzt[k] * w1[k / 2];
    
96.                 xzt[k + 1] = xzt[k + 1] * w1[k / 2];
    
97.                 sumXZTr += xzt[k]; sumXZTi += xzt[k + 1];
    
98.             }
    
99.             tempTr.Clear(); tempTi.Clear();
    
100.             for (int k = 0; k < 2 * iN; k += 2)
     
101.             {
     
102.                 tempTr.Add(xzt[k] - sumXZTr / iN);
     
103.                 tempTi.Add(xzt[k + 1] - sumXZTi / iN);
     
104.             }
     
105.             kbfft(tempTr.Count, 0, 0);
     
106.             for (int k = 0; k < iN / 2; k++)
     
107.             {
     
108.                 double temR = Convert.ToDouble(tempFr[k]), temI = Convert.ToDouble(tempFi[k]);
     
109.                 xz[k] += Math.Pow(Math.Sqrt(temI * temI + temR * temR) / iN * 2, 2);
     
110.             }
     
111.         }
     
112.         for (int i = 0; i < iN / 2; i++)
     
113.         {
     
114.             m_ZoomFFTSeries.Add(Math.Sqrt(xz[i] / iL));
     
115.         }
     
116.         m_ZoomXmin = fl; m_ZoomXmax = fh;
     
117.         m_Xstep = (fh - fl) / iN * 2;
     
118.     }
     
119.    
     
120.     /// <summary>
     
121.     /// GetZoomFFTSeries(double fin, int scale, int iN, int iL, int iM)的重构方法。
     
122.     /// </summary>
     
123.     public void GetZoomFFTSeries()
     
124.     {
     
125.         this.GetZoomFFTSeries(500, 5, 1024, 1, m_TimeSeries.Count/16);
     
126.     }
     
127.     #endregion
     
128. 
     
129.     #region 快速傅里叶变换方法
     
130.     /***************************************************************************
     
131.      * 入口参数：
     
132.      *  l: l = 0, 傅立叶变换; l = 1, 逆傅立叶变换
     
133.      * il: il = 0,不计算傅立叶变换或逆变换模和幅角；il = 1,计算模和幅角
     
134.      * n: 输入的点数，为偶数，一般为32，64，128，...,1024等
     
135.      * k: 满足n=2^k(k>0),实质上k是n个采样数据可以分解为偶次幂和奇次幂的次数
     
136.      * pr[]: l=0时，存放N点采样数据的实部
     
137.      * l=1时, 存放傅立叶变换的N个实部
     
138.      * pi[]: l=0时，存放N点采样数据的虚部
     
139.      * l=1时, 存放傅立叶变换的N个虚部
     
140.      *
     
141.      * 出口参数：
     
142.      * fr[]: l=0, 返回傅立叶变换的实部
     
143.      * l=1, 返回逆傅立叶变换的实部
     
144.      * fi[]: l=0, 返回傅立叶变换的虚部
     
145.      * l=1, 返回逆傅立叶变换的虚部
     
146.      * pr[]: il = 1,l = 0 时，返回傅立叶变换的模
     
147.      * il = 1,l = 1 时，返回逆傅立叶变换的模
     
148.      * pi[]: il = 1,l = 0 时，返回傅立叶变换的辐角
     
149.      * il = 1,l = 1 时，返回逆傅立叶变换的辐角
     
150.      ***************************************************************************/
     
151.     /// <summary>
     
152.     /// 快速傅里叶变换方法。
     
153.     /// </summary>
     
154.     void kbfft(int n,int l, int il)
     
155.     {
     
156.         int k = (int)Math.Log(n, 2);
     
157.         double[] pr = new double[n], pi = new double[n], fr = new double[n], fi = new double[n];
     
158.         for (int i = 0; i < n; i++)
     
159.         {
     
160.             if (i < tempTr.Count)
     
161.                 pr[i] = Convert.ToDouble(tempTr[i]);
     
162.             else
     
163.                 pr[i] = 0;
     
164.             if (i < tempTi.Count)
     
165.                 pi[i] = Convert.ToDouble(tempTi[i]);
     
166.             else
     
167.                 pi[i] = 0;
     
168.         }
     
169. 
     
170.         int  m, is1, j, nv;
     
171.         double p, q, s, vr, vi, poddr, poddi;
     
172. 
     
173.         //排序
     
174.         for (int it = 0; it <= n - 1; it++)
     
175.         {
     
176.             m = it; is1 = 0;
     
177.             for (int i = 0; i <= k - 1; i++)
     
178.             {
     
179.                 j = m / 2; is1 = 2 * is1 + (m - 2 * j); m = j;
     
180.                 fr[it] = pr[is1]; fi[it] = pi[is1];
     
181.             }
     
182.         }
     
183.         
     
184.         //蝶形运算
     
185.         pr[0] = 1.0; pi[0] = 0.0;
     
186.         p = Math.PI * 2 / (1.0 * n);
     
187.         pr[1] = Math.Cos(p); pi[1] = -Math.Sin(p);
     
188.         if (l != 0) pi[1] = -pi[1];
     
189.         for (int i = 2; i <= n - 1; i++)
     
190.         {
     
191.             p = pr[i - 1] * pr[1]; q = pi[i - 1] * pi[1];
     
192.             s = (pr[i - 1] + pi[i - 1]) * (pr[1] + pi[1]);
     
193.             pr[i] = p - q; pi[i] = s - p - q;
     
194.         }
     
195.         for (int it = 0; it <= n - 2; it = it + 2)
     
196.         {
     
197.             vr = fr[it]; vi = fi[it];
     
198.             fr[it] = vr + fr[it + 1]; fi[it] = vi + fi[it + 1];
     
199.             fr[it + 1] = vr - fr[it + 1]; fi[it + 1] = vi - fi[it + 1];
     
200.         }
     
201.         m = n / 2; nv = 2;
     
202.         for (int l0 = k - 2; l0 >= 0; l0--)
     
203.         {
     
204.             m = m / 2; nv = 2 * nv;
     
205.             for (int it = 0; it <= (m - 1) * nv; it = it + nv)
     
206.                 for ( j = 0; j <= (nv / 2) - 1; j++)
     
207.                 {
     
208.                     p = pr[m * j] * fr[it + j + nv / 2];
     
209.                     q = pi[m * j] * fi[it + j + nv / 2];
     
210.                     s = pr[m * j] + pi[m * j];
     
211.                     s = s * (fr[it + j + nv / 2] + fi[it + j + nv / 2]);
     
212.                     poddr = p - q; poddi = s - p - q;
     
213.                     fr[it + j + nv / 2] = fr[it + j] - poddr;
     
214.                     fi[it + j + nv / 2] = fi[it + j] - poddi;
     
215.                     fr[it + j] = fr[it + j] + poddr;
     
216.                     fi[it + j] = fi[it + j] + poddi;
     
217.                 }
     
218.         }
     
219.         if (l != 0)
     
220.             for (int i = 0; i <= n - 1; i++)
     
221.             {
     
222.                 fr[i] = fr[i] / (1.0 * n);
     
223.                 fi[i] = fi[i] / (1.0 * n);
     
224.             }
     
225.         if (il != 0)
     
226.             for (int i = 0; i <= n - 1; i++)
     
227.             {
     
228.                 pr[i] = Math.Sqrt(fr[i] * fr[i] + fi[i] * fi[i]);
     
229.                 pr[i] = (pr[i] / (n / 2)); //各次谐波幅值，其中pr[1]为基波幅值
     
230.                 if (Math.Abs(fr[i]) < 0.000001 * Math.Abs(fi[i]))
     
231.                 {
     
232.                     if ((fi[i] * fr[i]) > 0) pi[i] = 90.0;
     
233.                     else pi[i] = -90.0;
     
234.                 }
     
235.                 else
     
236.                     pi[i] = Math.Atan(fi[i] / fr[i]) * 360.0 / (2*Math.PI);
     
237.             }
     
238. 
     
239. 
     
240.         tempFi.Clear(); tempFr.Clear();
     
241.         for (int i = 0; i < n; i++)
     
242.         {
     
243.             tempFr.Add(fr[i]);
     
244.             tempFi.Add(fi[i]);
     
245.         }
     
246.     }
     
247.     #endregion

复制代码

这是我写的一个类（其中一部分）。。。希望对大家有所帮助