# 留言板  涂鸦板 lei5916418 2015-9-25 16:13 chunmu126 2012-12-25 11:31

1、改进的全相位时移相位差频谱分析算法
http://home.chinavib.com/blog-62061-18598.html

2、在下面的关于悬臂梁的分析中 chunmu126 2012-6-14 17:01

），所以我想对于基频的频率量先按照50hz来进行矫正，希望得到真正的基频频率大小，然后再进行最终的谐波计算，但是在这个过程中，若取ee = mod((p1-p2)/180/(1-1/N),Fs/N)，能够通过50hz矫正得到近似的真正的频率结果51.7491，但是若取ee = mod((p1-p2)/180/(1-1/N),1)，就不能得到正确的结果，具体如下：

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
…… chunmu126 2012-6-12 11:33

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 2012-1-12 11:41

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)

subplot(211);plot(a1.^2/N);
subplot(212);plot(a2);

2012-1-12 11:20 上传下载附件 (16.61 KB)

FFT的振幅与apFFT的振幅相近；但FFT的功率谱与apFFT的振幅根本就不在一个数量级上。

1.这是不是说明apFFT的振幅谱不能反映FFT的功率谱呢？还是说要在满足什么特定的条件下才能反映？
2.我举得例子里FFT的振幅谱怎么和apFFT功率谱一样呢，从数学原理上

这两个问题困扰了我很久，希望得到您的解答。谢谢 hjj73 2010-6-1 10:25  lsf_xaut 2009-5-13 13:48

我正在学习您与黄翔东老师的专著<数字信号全相位谱分析和滤波技术>一书，遇见几个自己不能解决的问题，现在列举如下，希望您能给与指导！
我的问题是关于课本中附录2和附录3的matlab程序问题.
关于附录2（FFT/apFFT综合校正程序），其中有语句
ee=mod((p1-p2)/180/(1-1/N),1);%频率偏离校正值

关于附录3(apFFT时移相位差法校正程序)，有语句
h=2*pi*g.*(1-g.*g)./sin(pi*g);
aa1=abs((h.^2).*a2)/2;

aa1=abs(a2./(h.^2));(根据式(4-58)得来的)
这两个问题我一直弄不明白，希望老师给予指导! lsf_xaut 2009-5-13 13:29

GMT+8, 2023-12-7 01:22 , Processed in 0.041555 second(s), 10 queries , Gzip On.