论坛 › 信号处理方法 › 加窗插值为何不准确，哪里出问题了？

dsfire 发表于 2010-1-25 11:06

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

这是我写的加窗插值代码，采用的是海明窗，修正公式如下图所示，
为什么测出来的幅值，频率和相位都差距那么大呢？
望赐教！
clc;
clear;
close all;
format short g;
f0=50.5;
fs=50*256;
N=2048;%2048/256=8
n=;
xb=;
Q=;
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;
for i=fix(35/(fs/N)):fix(65/(fs/N))
if u(i+1)>y
y=u(i+1)
num=i;
end
end
%★★★★★★谱线最大值已经找到，下面是要确定k1和k2还有y1和y2★★★★★★
k1=zeros(1,63);
k2=zeros(1,63);
k1(1)=num;
y1=u(k1(1)+1);%
k2(1)=k1(1)+1;
y2=u(k2(1)+1);

y3=u(num-1);
max=y2;
if y3>y2
max=y3;
end
if max==y3;
t=y1;
y1=max;
y2=t;
k1(1)=k1(1)-1;%根据公式，谱线谁在前谁是y1，对应的就是k1
k2(1)=k1(1)+1;%根据公式，谱线谁在后谁是y2，对应的就是k2
end
%%★★★★★★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)+a+0.5)*fs/N; %基波频率
Ph(1)=phase(v(k1(1)))+pi/2-pi*(a+0.5); %基波相位
%%★★★★★★到这里，基波频率，基波幅值，基波相位确定完毕★★★★★★


%%%%找出2到63次谐波的幅值，频率和相位
for i=2:63
k1(i)=fix(i*ff(1)/(fs/N));
k2(i)=k1(i)+1;
y1=u(k1(i)+1);
y2=u(k2(i)+1);
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; %各次谐波幅值，到63次谐波
ff(i)=(k1(i)+a+0.5)*fs/N; %各次谐波频率，到63次谐波
Ph(i)=phase(v(k1(i)))+pi/2-pi*(a+0.5); %各次谐波相位，到63次谐波
end
A
ff
Ph

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

鱼儿闯江湖 发表于 2010-1-27 17:12

回复 楼主 dsfire 的帖子

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

dsfire 发表于 2010-1-27 22:16

回复 沙发 鱼儿闯江湖 的帖子

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

songzy41 发表于 2010-1-30 07:03

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

dsfire 发表于 2010-1-30 19:48

回复 地板 songzy41 的帖子

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

songzy41 发表于 2010-1-31 09:08

我修改基波程序如下，谐波部分LZ自已再加上去：
clc;
clear;
close all;
format short g;
f0=50.5;
fs=50*256;
N=2048;%2048/256=8
n=;
xb=;
Q=;
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;
=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 发表于 2010-1-31 14:40

回复 6楼 songzy41 的帖子

您好，谢谢您的回复！
不知道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.368044.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.424780.96673   -5.616    1.0205   -4.7237    -5.1644
频率中偶次谐波测出来的差距比较大。
测出的谐波相角差距还是很大，我上面的程序哪里写错了么？
谢谢大侠了！

songzy41 发表于 2010-1-31 14:52

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

dsfire 发表于 2010-1-31 15:01

回复 8楼 songzy41 的帖子

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

songzy41 发表于 2010-1-31 18:29

原帖由 dsfire 于 2010-1-31 15:01 发表 http://www.chinavib.com/forum/images/common/back.gif
多谢回复，明白了！
那后面谐波频率和相位的差距还望您给指教！
对于计算的结果确实就是如此，结果并不好，说明用海明窗在谐波分析中并不理想。
海明窗主瓣和旁瓣衰减关系是-43dB/oct，即每倍频程衰减43dB。基波和2次谐波正好是一个倍频程，但基波和2次谐波相差70dB，所以计算2次谐波肯定受到基波的影响，计算就不正确了。同样偶次谐波幅值都比较小，奇次谐波比较大（实际情况就是这样），由于海明窗的衰减小，奇次谐波就影响了偶次谐波的参数计算。建议用更好的窗函数，例如Blackman-Harris窗函数。

freeplus 发表于 2010-2-1 11:15

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

freeplus 发表于 2010-2-1 11:17

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

dsfire 发表于 2010-2-1 13:40

回复 10楼 songzy41 的帖子

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

dsfire 发表于 2010-2-1 13:45

回复 12楼 freeplus 的帖子

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

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

Minelin 发表于 2014-12-1 10:57

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

查看完整版本: 加窗插值为何不准确，哪里出问题了？