求教 关于全相位 FFT/apFFT矫正程序的问题
下面是我在王兆华老师的FFT/apFFT矫正程序的基础上，写的关于模拟电力谐波的分析程序，由于我信号处理方面的知识很少，所以还有很多不理解的地方，请王老师和各位高手帮我解决一下。close all;clc;clear all;
N=1024;%现在固定在个数据点数
f11=51.1;
Fs=2500;
n=-N+1:N-1;
apl=;
phol=;
%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)+; %构成长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();xlim(); grid
subplot(5,1,2);stem(a2,'.');title('apFFT amplitude spectrum');ylim();xlim(); grid
subplot(5,1,3);stem(p2,'.');title('apFFT phase spectrum');ylim();xlim(); grid
subplot(5,1,4);stem(ee,'.');title('frequency correction spectrum');ylim([-1,1]);xlim(); grid
subplot(5,1,5);stem(aa,'.');title('amplitude correction spectrum');ylim();xlim(); 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.0000305.0000350.5000 33.0000346.0000 56.0000 7.0000326.0000300.0000
Columns 10 through 18
285.0000296.0000282.0000280.0000285.0000295.0000290.0000343.0000355.0000
Columns 19 through 21
70.0000290.0000312.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，不知道该怎么改正，请王老、各位高手帮我解答一下，谢谢！
本帖最后由 zhwang554 于 2012-6-12 22:10 编辑
回复 1 # chunmu126 的帖子
a按照书中P62页的公式（4-45）得到的是角频率, 除以2*pi/N得到的频率
对1进行取余数是使计算频偏 >0, 以适应下面计算频率校正值公式
b Matlab中数椐从1开始,计算得到的频谱频率从直流0开始
c将信号sin 改为cos即可, fft中的相位是对cos的
谢谢王老的的解答，明白了！ 本帖最后由 chunmu126 于 2012-6-14 17:24 编辑
[*]王兆华老师，谢谢您的答疑，我现在对于问题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_bodisp('频率校正值')(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_bodisp('频率校正值')(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取余数后在求绝对值，这三种之间有什么不同呢？烦请王老师帮我解答一下，谢谢！
本帖最后由 zhwang554 于 2012-6-15 02:44 编辑
回复 4 # chunmu126 的帖子
附录2和附录3的2个程序都是用於理解fft/apfft和 apfft/apfft校正原理的程序, 不是实用程序, 所以
1)采样频率选为N, 即fft的频率分辨率为1Hz/格, 不作频率转换使初学者搞不清楚
2)不用相位差判断语句, 所以出现振幅校正公式中有floor语句
3)不加噪音
4)不加自动搜索振幅峰值语句, rr=round(f1)是为了找出谱峰线的位置（rr+1），
5)加已知余弦信号, 这样可以知道校正是否正确
采样频率选为fs, 加自动搜索振幅峰值和加相位差判断语句的 apfft/apfft 和 fft/apfft 校正程序,
见本论坛zhwang554日志
”自动搜索峰值的 apfft/apfft 和 fft/apfft 校正程序”
http://forum.vibunion.com/home-spa ... -blog-id-18060.html
这个程序可分析实际数据，也可对已知频率(包恬振幅,相位)做估计
rr=round(i*f11)是为了找出谱峰线的位置，你设置rr=round(i*50), 这样对许多频率成份找出谱峰线的位置就不对了,笫1个基频就不对,实际峰值在rr=round(50.749123)=51,你设置成rr=round(1*50)=50, 频率校正值当然不对了. 不知(基频)频率值, 加自动搜索振幅峰值语句, 找出谱峰线的位置. 不能任意设置rr=round(i*50)你提出对频率分辨率来取模, xie_bo2参与时(这设置就不对)虽对笫一个频率是对的, 其它谐波都不对, 没有xie_bo2参与(本贴一楼的程序)得到的频率都不对,这有什么研究意义
所以找出谱峰线,对1取模用相位差判断语句时不取模对1进行取余数是使计算频偏在 0-1 之间 采样点数与采样频率有什么关系么？只采一个周期的是不是误差大？ jinyi7016 发表于 2014-10-23 16:08
采样点数与采样频率有什么关系么？只采一个周期的是不是误差大？
刚好采一个周期的话就是整周期采样，不用校正了。
一般来讲应该多采几个周期，以减小负频率的干涉。 整周期采样,但不是同步，仍然在泄漏吧。
我采样8周期，按以上方法，波形的初相位对计算结果影响挻大的，是什么原因？
页:
[1]