加窗插值为何不准确，哪里出问题了？

dsfire

这是我写的加窗插值代码，采用的是海明窗，修正公式如下图所示，
为什么测出来的幅值，频率和相位都差距那么大呢？
望赐教！

1. clc;
   
2. clear;
   
3. close all;
   
4. format short g;
   
5. f0=50.5;
   
6. fs=50*256;
   
7. N=2048;%2048/256=8
   
8. n=[0:1:N-1];
   
9. xb=[220,0.5,25,0.4,6,0.3,4,0.2,2];
   
10. Q=[pi/3,pi/3,pi/3,pi/3,pi/3,pi/3,pi/3,pi/3,pi/3];
    
11. s=zeros(1,N);
    
12. for i=1:9
    
13. s=s+xb(i)*cos(2*pi*f0*i*n./fs+Q(i));
    
14. end
    
15. figure (1);
    
16. plot(s);
    
17. 
    
18. w=0.54-0.46*cos(2*pi*n./N);
    
19. r=s.*w;
    
20. v=fft(r,N);
    
21. u=abs(v);
    
22. 
    
23. %为了求出谐波，那么必须找出基波，在35-65中频谱最大值所对应的就是离基波频率最近的那个点
    
24. y=0;
    
25. for i=fix(35/(fs/N)):fix(65/(fs/N))
    
26. if u(i+1)>y
    
27. y=u(i+1)
    
28. num=i;
    
29. end
    
30. end
    
31. %★★★★★★谱线最大值已经找到，下面是要确定k1和k2还有y1和y2★★★★★★
    
32. k1=zeros(1,63);
    
33. k2=zeros(1,63);
    
34. k1(1)=num;
    
35. y1=u(k1(1)+1);%
    
36. k2(1)=k1(1)+1;
    
37. y2=u(k2(1)+1);
    
38. 
    
39. y3=u(num-1);
    
40. max=y2;
    
41. if y3>y2
    
42. max=y3;
    
43. end
    
44. if max==y3;
    
45. t=y1;
    
46. y1=max;
    
47. y2=t;
    
48. k1(1)=k1(1)-1;%根据公式，谱线谁在前谁是y1，对应的就是k1
    
49. k2(1)=k1(1)+1;%根据公式，谱线谁在后谁是y2，对应的就是k2
    
50. end
    
51. %%★★★★★★k1和k2还有y1和y2确定完毕★★★★★★
    
52. 
    
53. %%★★★★★★确定αβ，在这里用a和b表示，确定基波频率★★★★★★
    
54. A=zeros(1,63);
    
55. ff=zeros(1,63);
    
56. Ph=zeros(1,63);
    
57. b=(y2-y1)/(y2+y1);
    
58. a=1.21874943*b+0.13349531*b^3+0.05301420*b^5+0.03656014*b^7;
    
59. A(1)=(y1+y2)*(2.26557103+1.22719978*a^2+0.37607775*a^4+0.09767389*a^6)/N; %基波幅值
    
60. ff(1)=(k1(1)+a+0.5)*fs/N; %基波频率
    
61. Ph(1)=phase(v(k1(1)))+pi/2-pi*(a+0.5); %基波相位
    
62. %%★★★★★★到这里，基波频率，基波幅值，基波相位确定完毕★★★★★★
    
63. 
    
64. 
    
65. %%%%找出2到63次谐波的幅值，频率和相位
    
66. for i=2:63
    
67. k1(i)=fix(i*ff(1)/(fs/N));
    
68. k2(i)=k1(i)+1;
    
69. y1=u(k1(i)+1);
    
70. y2=u(k2(i)+1);
    
71. b=(y2-y1)/(y2+y1);
    
72. a=1.21874943*b+0.13349531*b^3+0.05301420*b^5+0.03656014*b^7;
    
73. A(i)=(y1+y2)*(2.26557103+1.22719978*a^2+0.37607775*a^4+0.09767389*a^6)/N; %各次谐波幅值，到63次谐波
    
74. ff(i)=(k1(i)+a+0.5)*fs/N; %各次谐波频率，到63次谐波
    
75. Ph(i)=phase(v(k1(i)))+pi/2-pi*(a+0.5); %各次谐波相位，到63次谐波
    
76. end
    
77. A
    
78. ff
    
79. Ph

复制代码

[ 本帖最后由 dsfire 于 2010-1-25 11:08 编辑 ]

鱼儿闯江湖

:loveliness: 你好，我也做过你这方面的仿真，我的效果也没有论文上说的那么好，但是我觉得误差不是很大，频率的还是很大的，想讨论下，你的误差是多少的啊

dsfire

上面的程序一运行就知道了。相位误差比较大，我觉得是我写的相位校正的matlab式子有问题。频率我就是用上面的式子，误差也不小。不过，文章上说，先求出基波频率F，然后直接用i*F来代替各次谐波频率（i=2，3，4 ......）。我上面的程序不是这样做的。

songzy41

LZ是按哪一篇文献来编写的？

dsfire

应用FFT进行电力系统谐波分析的改进算法 --中国电机工程学报
其实有好几篇文章都有一样的校正式子。
望高手给予指点！

songzy41

我修改基波程序如下，谐波部分LZ自已再加上去：
clc;
clear;
close all;
format short g;
f0=50.5;
fs=50*256;
N=2048;%2048/256=8
n=[0:1:N-1];
xb=[220,0.5,25,0.4,6,0.3,4,0.2,2];
Q=[pi/3,pi/3,pi/3,pi/3,pi/3,pi/3,pi/3,pi/3,pi/3];
s=zeros(1,N);
for i=1:9
s=s+xb(i)*cos(2*pi*f0*i*n./fs+Q(i));
end
%figure (1);
%plot(s);

w=0.54-0.46*cos(2*pi*n./N);
r=s.*w;
v=fft(r,N);
u=abs(v);

%为了求出谐波，那么必须找出基波，在35-65中频谱最大值所对应的就是离基波频率最近的那个点
y=0;
n1=fix(35/(fs/N))+1; n2=fix(65/(fs/N))+1;
[y,num]=max(u(n1:n2));
num=num+n1-1;
%★★★★★★谱线最大值已经找到，下面是要确定k1和k2还有y1和y2★★★★★★
k1=zeros(1,63);
k2=zeros(1,63);
if u(num+1)<u(num-1)
    num=num-1;
end
k1(1)=num; k2(1)=num+1;
y1=u(num); y2=u(num+1);
%%★★★★★★k1和k2还有y1和y2确定完毕★★★★★★

%%★★★★★★确定αβ，在这里用a和b表示，确定基波频率★★★★★★
A=zeros(1,63);
ff=zeros(1,63);
Ph=zeros(1,63);
b=(y2-y1)/(y2+y1);
a=1.21874943*b+0.13349531*b^3+0.05301420*b^5+0.03656014*b^7;
A(1)=(y1+y2)*(2.26557103+1.22719978*a^2+0.37607775*a^4+0.09767389*a^6)/N; %基波幅值
ff(1)=(k1(1)-1+a+0.5)*fs/N; %基波频率
Ph(1)=phase(v(k1(1)))-pi*(a+0.5); %基波相位
%%★★★★★★到这里，基波频率，基波幅值，基波相位确定完毕★★★★★★
fprintf('%5.6f   %5.6f   %5.6f\n',A(1),ff(1),Ph(1));

dsfire

您好，谢谢您的回复！
不知道Ph(1)=phase(v(k1(1)))-pi*(a+0.5); %基波相位 这句中的pi/2为什么去掉呢?

下面是我谐波的程序
for i=2:9
    k1(i)=fix(i*ff(1)/(fs/N))+1;
    k2(i)=k1(i)+1;
    y1=u(k1(i));
    y2=u(k2(i));
    b=(y2-y1)/(y2+y1);
    a=1.21874943*b+0.13349531*b^3+0.05301420*b^5+0.03656014*b^7;
    A(i)=(y1+y2)*(2.26557103+1.22719978*a^2+0.37607775*a^4+0.09767389*a^6)/N; %各次谐波幅值，到9次谐波
    ff(i)=(k1(i)-1+a+0.5)*fs/N; %各次谐波频率，到9次谐波
    Ph(i)=phase(v(k1(i)))-pi*(a+0.5); %各次谐波相位，到9次谐波
end  
A(1:9)
ff(1:9)
Ph(1:9)
运行之后出来的数据：
次数： 基波           2           3                  4               5                  6             7                8              9           
幅值：219.91   0.56293   24.9          0.48755   5.9803    0.36804  4.0086   0.26062   2.0066
频率：50.505   106.95     151.58     206.32     252.71     305.25    353.64   404.73     454.53
相角：1.0443    -1.5182   1.0089    -0.42478  0.96673   -5.616    1.0205   -4.7237    -5.1644
频率中偶次谐波测出来的差距比较大。
测出的谐波相角差距还是很大，我上面的程序哪里写错了么？
谢谢大侠了！

songzy41

在LZ提供的文献“应用FFT进行电力系统谐波分析的改进算法”（中国电机工程学报）中是对正弦信号处理的（见式(1)），而对余弦信号（LZ提供的是余弦信号）差pi/2。

dsfire

多谢回复，明白了！
那后面谐波频率和相位的差距还望您给指教！

songzy41

原帖由 dsfire 于 2010-1-31 15:01 发表

                               
登录/注册后可看大图

多谢回复，明白了！
那后面谐波频率和相位的差距还望您给指教！

对于计算的结果确实就是如此，结果并不好，说明用海明窗在谐波分析中并不理想。
海明窗主瓣和旁瓣衰减关系是-43dB/oct，即每倍频程衰减43dB。基波和2次谐波正好是一个倍频程，但基波和2次谐波相差70dB，所以计算2次谐波肯定受到基波的影响，计算就不正确了。同样偶次谐波幅值都比较小，奇次谐波比较大（实际情况就是这样），由于海明窗的衰减小，奇次谐波就影响了偶次谐波的参数计算。建议用更好的窗函数，例如Blackman-Harris窗函数。

freeplus

呵呵，楼主，我没有看过论文，粗粗看了一下你的源代码，发现在确定y1,y2,y3幅值比较中有一条明显错误语句：
k1(1)=k1(1)-1;%根据公式，谱线谁在前谁是y1，对应的就是k1
k2(1)=k1(1)+1;%根据公式，谱线谁在后谁是y2，对应的就是k2
           ^^^^^
是不是 k2(1)=k2(1)+1; ???

freeplus

另外说一句，你的代码实在太C语言化了，很多循环语句都可以用MATLAB并行特性简化为一条语句的。

dsfire

非常感谢您的回复，我再仔细的看看！

dsfire

哈哈，我刚学matlab语言，对里面的语句还没有太多认识！
不过，那两个语句是应该是对的，因为k1和k2是顺序的，k2必须等于k1+1。
同样，非常感谢你的回复！
以后有什么问题也希望你能指教！

[ 本帖最后由 dsfire 于 2010-2-1 14:47 编辑 ]

Minelin

学习了，刚学有关FFT这方面的知识。。。