求教 关于全相位 FFT/apFFT矫正程序的问题

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，不知道该怎么改正，请王老、各位高手帮我解答一下，谢谢！
  

zhwang554

回复 1 # chunmu126 的帖子

a  按照书中P62页的公式（4-45）得到的是角频率, 除以2*pi/N得到的频率
                             
     对1进行取余数是使计算频偏 >0, 以适应下面计算频率校正值公式

b   Matlab中数椐从1开始,计算得到的频谱频率从直流0开始

c  将信号sin 改为cos即可, fft中的相位是对cos的

chunmu126

谢谢王老的的解答，明白了！

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取余数后在求绝对值，这三种之间有什么不同呢？烦请王老师帮我解答一下，谢谢！
  
  

zhwang554

回复 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

采样点数与采样频率有什么关系么？只采一个周期的是不是误差大？

yangzj

jinyi7016 发表于 2014-10-23 16:08
采样点数与采样频率有什么关系么？只采一个周期的是不是误差大？

刚好采一个周期的话就是整周期采样，不用校正了。
一般来讲应该多采几个周期，以减小负频率的干涉。

jinyi7016

整周期采样,但不是同步，仍然在泄漏吧。
我采样8周期，按以上方法，波形的初相位对计算结果影响挻大的，是什么原因？