留言板
- lei5916418
- 王老师你好，我最近看了您的apfft程序和相关文献，有几个问题想向您请教
使用hanning卷积窗的振幅校正公式为：AA(m)=2*((1-dk.*dk)./sinc(dk))^2.*abs(peak_F1(m));这个我能理解
按照我的推导，使用矩形卷积窗的振幅校正公式应当为AA(m)=2/(2*sinc(dk))^2.*abs(peak_F1(m));
可是实际计算中AA(m)=2/(sinc(dk))^2.*abs(peak_F1(m));才是正确的，两者之间差了整4倍，这个是我不能理解的；
另外，如果使用hamming卷积窗，那么校正公式应当为AA(m)=2*((1-dk.*dk)./sinc(dk)./(1.08-0.16*dk^2))^2.*abs(peak_F1(m));
可是实际计算结果完全不正确，而我也没有找到正确的hamming卷积窗的振幅校正结果，所以我想向您请教apfft下振幅校正公式的推导过程，能否请您以hamming窗为例向我讲解一下呢
另外说明，我所使用的程序来自http://home.vibunion.com/blog-62061-18060.html这篇文章
- chunmu126
- 王兆华老师，您提出的全相位算法的基础是，可以准确的测量信号的初始相位。经过对王兆华老师及网友分享的代码的实验发现，如果要对实际采样的得到的信号进行直接运算的话，不能得到其真实的初始相位，但复制和频率的校正结果还是很好的。
例如：实际得到的信号，不可能从t=0开始，也不会有t<0的部分，王老师在教学视频中也提到了这一点。但是在实际的仿真中，我分别使用fft/apfft法、apfft/apfft法，以"t=1:2*N-1"、"t=1:3*N-1"来合成余弦信号序列进行仿真。但是结果，除了可以获得高精度的幅值和频率校正结果外（与选择"t=-N+1:N-1"、"t=-N+1:2*N-1"来合成余弦信号序列的精度基本相同），不能得到它的真实相位，而且还差的很远。
我知道这里之所以得到的幅值和频率的校正结果，精度很高，是因为全相位时移相位差法在这种情况下是可用的。也就是王老师所说的，全相位和截取序列的起点是没有关系的。
但是，我现在不知道如何才能“校正”出真实的相位来。需要对代码如何改进？
请王老师和各位网友给我指导一下，谢谢！
补充：
1、改进的全相位时移相位差频谱分析算法
http://home.chinavib.com/blog-62061-18598.html
文中提到了对选取中心样本点的改进，但是文中也提到了，是以牺牲精度为代价的，所以王老师一定有更好的办法。
2、在下面的关于悬臂梁的分析中
http://www.chinavib.com/forum.php?mod=viewthread&tid=74297&page=1#pid416482
我也没有能够得到如王老师，算出的相位值来，而且差别很大，请王老师，给与指导。
- chunmu126
- 王兆华老师，谢谢您的答疑，我现在对于问题a、中为在求频率偏移量的时候选择用1，来求余数，我还是有疑问，请您帮我解答一下，谢谢！
事情是这样的，因为在实际当中，我不知道采集到的信号它的基频是多少（即程序中的f11=51.749123
），所以我想对于基频的频率量先按照50hz来进行矫正，希望得到真正的基频频率大小，然后再进行最终的谐波计算，但是在这个过程中，若取ee = mod((p1-p2)/180/(1-1/N),Fs/N)，能够通过50hz矫正得到近似的真正的频率结果51.7491，但是若取ee = mod((p1-p2)/180/(1-1/N),1)，就不能得到正确的结果，具体如下：
程序的主体与一楼的一样。
首先声明：xie_bo=zeros(1,21);xie_bo2=zeros(1,21);
for i=1:21
y=y+apl(i)*cos(2*pi*i*f11*n/Fs+phol(i)*pi/180);
xie_bo(i)=i*f11; %模拟采集到的各次谐波的频率
xie_bo2(i)=i*50;%按照50hz得到的各次谐波频率
end
若采用
ee = mod((p1-p2)/180/(1-1/N),1);
r=round(xie_bo2*N/Fs)
xie_bo
disp('频率校正值')
(ee(r+1)+floor(xie_bo2*N/Fs))*Fs/N
则无法得到较为准确的基频，结果如下：
xie_bo =
Columns 1 through 4
51.749123
103.498246
155.247369
206.996492
Columns 5 through 8
258.745615
310.494738
362.243861
413.992984
Columns 9 through 12
……
频率校正值
ans =
Columns 1 through 4
49.3077113000786
98.6186179056816
150.359933774583
199.67687346012
Columns 5 through 8
251.413025795446
298.265414236582
350.024879021602
399.306840999873
Columns 9 through 12
……
但是若采用
ee = mod((p1-p2)/180/(1-1/N),Fs/N);
r=round(xie_bo2*N/Fs)
xie_bo
disp('频率校正值')
(ee(r+1)+floor(xie_bo2*N/Fs))*Fs/N
则能够得到较为准确的基频，结果如下：
xie_bo =
Columns 1 through 4
51.749123
103.498246
155.247369
206.996492
Columns 5 through 8
258.745615
310.494738
362.243861
413.992984
Columns 9 through 12
……
频率校正值
ans =
Columns 1 through 4
51.7491175500786
101.060024155682
153.878992002122
203.195931687659
Columns 5 through 8
254.932084022985
300.706820486582
352.466285271602
401.748247249873
Columns 9 through 12
……
从实验的结果看，我就无法理解了为什么不是对频率分辨率来取模了，
但是对于本贴一楼的程序来说（即没有xie_bo2的参与），如果选择对频率分辨率Fs/N来取模的话，得到的频率却不对了，所以现在很迷茫，请王老师帮我解答一下！
还有就是如果“对1进行取余数是使计算频偏 >0, 以适应下面计算频率校正值公式”，那么对其他的正数取余数，或对Fs/N取余数后在求绝对值，这三种之间有什么不同呢？烦请王老师帮我解答一下，谢谢！
- chunmu126
- 下面是我在王兆华老师的FFT/apFFT矫正程序的基础上，写的关于模拟电力谐波的分析程序，由于我信号处理方面的知识很少，所以还有很多不理解的地方，请王老师帮我解决一下。
close all;clc;clear all;
N=1024;%现在固定在个数据点数
f11=51.1;
Fs=2500;
n=-N+1:N-1;
apl=[220.5,4.4,10,3,6,1,3,0.8,1.5,0.8,1.1,0.04,0.85,0.1,1,0.04,0.5,0.02,0.3,0.005,0.01];
phol=[10,35,80.5,123,76,146,97,56,30,15,26,12,10,15,25,20,73,85,160,20,42];
%phol=zeros(1,21);%如果相位全选为0，则相位矫正结果全为270
ci=21;
y=zeros(1,2*N-1);
xie_bo=zeros(1,21);
for i=1:ci
y=y+apl(i)*sin(2*pi*i*f11*n/Fs+phol(i)*pi/180);%各次谐波叠加 i*f11
xie_bo(i)=i*f11;%计算谐波频率
end
%先进行传统的FFT
y1 = y(N:2*N-1);
win = hanning(N)';
win1 = win/sum(win); %窗归一 用于把加权求和的窗进行归一化
y11 = y1.*win1;
y11_fft = fft(y11,N);%做N点的FFT运算
a1 = abs(y11_fft); %FFT振幅谱
p1 = mod(phase(y11_fft)*180/pi,360); %FFT相位谱 angle(v)也可以求
%再进行apFFT进行运算
y2 = y(1:2*N-1); %2N-1个输入数据
win = hanning(N)';
winn = conv(win,win); %apFFT须要卷积窗
win2 = winn/sum(winn); %窗归一
y22 = y2.*win2;
y222 = y22(N:end)+[0 y22(1:N-1)]; %构成长N的apFFT输入数据 求和
y2_fft = fft(y222,N);
a2 = abs(y2_fft); %apFFT振幅谱
p2 = mod(phase(y2_fft)*180/pi,360);
%校正量的计算
ee = mod((p1-p2)/180/(1-1/N),1); %频率偏离矫正值
aa = (a1.^2)./a2*2; %振幅校正值 为什么要乘以2
subplot(5,1,1);stem(a1,'.');title('FFT amplitude spectrum');ylim([0,1]);xlim([0 N/2]); grid
subplot(5,1,2);stem(a2,'.');title('apFFT amplitude spectrum');ylim([0,110]);xlim([0 N/2]); grid
subplot(5,1,3);stem(p2,'.');title('apFFT phase spectrum');ylim([0,400]);xlim([0 N/2]); grid
subplot(5,1,4);stem(ee,'.');title('frequency correction spectrum');ylim([-1,1]);xlim([0 N/2]); grid
subplot(5,1,5);stem(aa,'.');title('amplitude correction spectrum');ylim([0,420]);xlim([0 N/2]); grid
xie_bo
r=round(xie_bo*N/Fs)%求谱线个数
disp('相位校正值')
p2(r+1)
disp('频率校正值')
(ee(r+1)+floor(xie_bo*N/Fs))*Fs/N %
disp('振幅校正值')
aa(r+1)
得到的结果如下：
xie_bo =
1.0e+003 *
Columns 1 through 9
0.0511 0.1022 0.1533 0.2044 0.2555 0.3066 0.3577 0.4088 0.4599
Columns 10 through 18
0.5110 0.5621 0.6132 0.6643 0.7154 0.7665 0.8176 0.8687 0.9198
Columns 19 through 21
0.9709 1.0220 1.0731
r =
Columns 1 through 16
21 42 63 84 105 126 147 167 188 209 230 251 272 293 314 335
Columns 17 through 21
356 377 398 419 440
相位校正值
ans =
Columns 1 through 9
280.0000 305.0000 350.5000 33.0000 346.0000 56.0000 7.0000 326.0000 300.0000
Columns 10 through 18
285.0000 296.0000 282.0000 280.0000 285.0000 295.0000 290.0000 343.0000 355.0000
Columns 19 through 21
70.0000 290.0000 312.0000
频率校正值
ans =
1.0e+003 *
Columns 1 through 9
0.0511 0.1022 0.1533 0.2044 0.2555 0.3066 0.3577 0.4088 0.4599
Columns 10 through 18
0.5110 0.5621 0.6132 0.6643 0.7154 0.7665 0.8176 0.8687 0.9198
Columns 19 through 21
0.9709 1.0220 1.0731
振幅校正值
ans =
Columns 1 through 9
220.5000 4.3962 9.9995 2.9997 5.9998 0.9998 2.9999 0.8001 1.5000
Columns 10 through 18
0.8000 1.1000 0.0400 0.8500 0.1000 1.0000 0.0400 0.5000 0.0200
Columns 19 through 21
0.3000 0.0050 0.0100
问题：
a、程序中关于频率偏移量的计算ee = mod((p1-p2)/180/(1-1/N),1); 应该是按照书中P62页的公式（4-45）得到的，但是我从公式中推导的程序为ee=mod( ( (p1-p2)*pi/180 )*2/(n-1),1 );不知道这个对不对，王老师的表达式又是怎么得来的呢？另外，为什么是对1进行取余数呢，这个与频率分辨率有关系么？请王老帮我解答一下，谢谢！
b、程序中关于矫正结果的提取都用到了（r+1），这个是为什么呢，不应该是第r跟谱线的频率值，加上矫正量么？请王老帮我解答一下，谢谢！
c、从我的程序的结果可以看出，相位的矫正结果不对，我的不明白是为什么。如果把21个初始的相位全都定为0，则可到的相位矫正结果全是270，不知道该怎么改正，请王老帮我解答一下，谢谢！
- Oboyer
- 老师你好：
资料上说：apFFT的振幅谱就反映了FFT的功率谱，我用您书里附录3的“fft/apfft综合校正法”程序做了验证，发现一些问题，希望得到您的指教，程序如下（我只取了程序里的信号与求振幅的部分）:
N=128;
t=-N+1:N-1;
f1=19.333333333333333;A1=0.888888888888888;ph1=100.111111111111111111;
y=A1*exp(j*(2*pi*t*f1/N+ph1*pi/180));
y1 = y(N:2*N-1);%后N个输入数据
win = hanning(N)';;
win1 = win/sum(win);%窗归1
y11= y1.*win1;
y11_fft = fft(y11,N);
a1 = abs(y11_fft);%FFT振幅谱
y2 = y(1:2*N-1);%2N-1个输入数据
winn = conv(win,win);%apFFT须要卷积窗
win2 = winn/sum(winn);%窗归1
y22= y2.*win2;
y222=y22(N:end)+[0 y22(1:N-1)];%构成长N的apFFT输入数据
y2_fft = fft(y222,N);;
a2 = abs(y2_fft);%apFFT振幅谱
subplot(211);plot(a1);
subplot(212);plot(a2);
2012-1-12 11:15 上传下载附件 (17.14 KB)
第一幅图是FFT的振幅；第二幅图是apFFT的振幅。
下面绘制FFT的功率谱与apFFT的振幅谱：
subplot(211);plot(a1.^2/N);
subplot(212);plot(a2);
2012-1-12 11:20 上传下载附件 (16.61 KB)
第一幅图是FFT的功率谱；第二幅图是apFFT的振幅。
由以上四幅图明显可以看出：
FFT的振幅与apFFT的振幅相近；但FFT的功率谱与apFFT的振幅根本就不在一个数量级上。
我的问题是：
1.这是不是说明apFFT的振幅谱不能反映FFT的功率谱呢？还是说要在满足什么特定的条件下才能反映？
2.我举得例子里FFT的振幅谱怎么和apFFT功率谱一样呢，从数学原理上
推导，apFFT的振幅不应该是FFT振幅的平方才对么？
这两个问题困扰了我很久，希望得到您的解答。谢谢
- hjj73
- 你好，我在电厂从事设备管理，您是教师吗？
- 李元
- 王老师好
- lsf_xaut
- 王老师，您好！
我正在学习您与黄翔东老师的专著<数字信号全相位谱分析和滤波技术>一书，遇见几个自己不能解决的问题，现在列举如下，希望您能给与指导！
我的问题是关于课本中附录2和附录3的matlab程序问题.
关于附录2（FFT/apFFT综合校正程序），其中有语句
ee=mod((p1-p2)/180/(1-1/N),1);%频率偏离校正值
对于这个校正公式，我提出这样的问题：该语句应该为ee=mod((p1-p2)/180/(N-1),1);(这是按照课本P62式（4-45）得到的，它应该除以N-1),而在附录2中却是（1-1/N），这一点我很不明白！
关于附录3(apFFT时移相位差法校正程序)，有语句
h=2*pi*g.*(1-g.*g)./sin(pi*g);
aa1=abs((h.^2).*a2)/2;
要是按照课本（4-8）（4-43）（4-58），我个人觉得应该将该语句改为如下形式：
aa1=abs(a2./(h.^2));(根据式(4-58)得来的)
这两个问题我一直弄不明白，希望老师给予指导!
- lsf_xaut
- 您好！
