[已解决] 全相位FFT相位不正确

[复制链接]
Ljgupup 发表于 2021-6-25 19:55:28
本帖最后由 Ljgupup 于 2021-6-25 19:59 编辑

我最近学习王兆华老师的全相位FFT(apfft)算法,此算法需要输入2N-1个数据(本例N=256),测试发现只有当采样间隔t=(-nfft+1):(nfft-1);时画出来的相位谱显示的相位值和理论的值一样(本例是100°),但是我把采样间隔设置为t=0:(2*N-2);(总计2N-1个点)时,画出来的相位谱对应谱线的相位值完全不对。很奇怪,除了采样间隔设置不一样,其他地方都一样,搞不懂为什么,看王兆华老师的输入值的定义域为n∈(-N+1, N-1),我就想不通为什么要这样定义。

哎,真是想不明白呀,图二正确说明算法正确呀,而且t=0:(2*nfft-2)这样采样我是想模拟实际采样嘛,实际分析时数据就是连续的2*N-1个离散点呀?求大佬指点呀!!

%全相傅里叶变换
clear
clc
close all hidden
%准备数据
sf=256; %采样频率
nfft=256; %FFT计算点数
n=2*nfft-1; %数据采样点数

%采样点
t=(-nfft+1):(nfft-1);
% t=0:1:(2*nfft-2);
% t=1:(2*nfft-1);

x=0.5*cos(2*pi*17.92*t/sf+100*pi/180); %采样函数频率17.92,相位100°

% x=cos(50.3*2*pi*t/fs+pi/6)+cos(85.3*2*pi*t/sf+pi/3)+cos(121.4*2*pi*t/sf+pi/2);
% x=0.5*cos(2*pi*50*t/sf+pi/6);


% =============全相位FFT预处理开始(双窗)================
win1=hanning(nfft)';
win2=hanning(nfft)';
win=conv(win1,win2); %卷积窗
win=win/sum(win);     %归一化
w=pi*2;
y1=x.*win;                %加权
y1a=y1(nfft:end) + [0 y1(1:nfft-1)];
%=============全相位FFT预处理结束================

%计算FFT
Y=fft(y1a, nfft); %apFFT
Y1=fft(x(1:nfft),nfft); %传统FFT

subplot(2,2,1);
a1=abs(Y1)/nfft*2;
a1(1)=a1(1)*2;
plot(a1); %幅度谱
title('传统FFT');
subplot(2,2,2);
a=sqrt(abs(Y));
plot(a); %幅度谱
title('apFFT');
subplot(2,2,3);
plot(angle(Y1)*180/pi); %相位谱
subplot(2,2,4);
plot(angle(Y)*180/pi); %相位谱

复制代码



下载 (1).png
下载.png

最佳答案


zhwang554 发表于 2021-7-13 03:42:06
本帖最后由 zhwang554 于 2021-7-13 21:39 编辑

FFT相位谱中的相位值是序到(t=0:N-1)第一个采样点(t=0)的余弦相位值,初相位x=0.5*cos(2*pi*17.92*t/sf+100*pi/180)中的那个100。

apfft序列共有2N-1个样点,经予处理后形成N个样点再作FFT,予处理后N个样点的第一个样点(t=0)是2N-1个序列中间的那个样点,所以apfft的相位谱显示的是2N-1序列[t=(-nfft+1):(nfft-1)]中间点(t=0)的相位,而不是2N-1序列[t=0:(2*N-2)]第一点的相位。

补充内容 (2021-11-20 22:01):
楼主认为图2显示的相位46度不正确。图2是t=0:2N-2吋的相位图,它显示的不是初相100度,而是序列第N点的相位46度。可以计算出当初相100时的第N点相46度。

补充内容 (2021-11-20 22:09):
楼主认为图2的相位图显示的46度不正确。图2是t=0:2*N-2时的apfft相位图, 它显示这个序列第N点的相位值是46度。这个序列的第一个样点是100度,笫N点样点的相位就是46度,没有錯。

补充内容 (2021-11-20 22:11):
没有錯

补充内容 (2021-11-21 21:35):
为什么”这个序列的第一个样点相位是100度,第N个相位是46度”
請阅读下面的贴子。有详细的计算过程。

21 条回复


Ljgupup 发表于 2021-6-26 08:07:40
没人看么

Ljgupup 发表于 2021-6-26 16:37:23
本帖最后由 Ljgupup 于 2021-6-28 08:21 编辑

顶一个!!!!!!

Ljgupup 发表于 2021-6-30 14:53:27
顶一个!!!!!!

zhwang554 发表于 2021-7-13 03:42:06
本帖最后由 zhwang554 于 2021-7-13 21:39 编辑

FFT相位谱中的相位值是序到(t=0:N-1)第一个采样点(t=0)的余弦相位值,初相位x=0.5*cos(2*pi*17.92*t/sf+100*pi/180)中的那个100。

apfft序列共有2N-1个样点,经予处理后形成N个样点再作FFT,予处理后N个样点的第一个样点(t=0)是2N-1个序列中间的那个样点,所以apfft的相位谱显示的是2N-1序列[t=(-nfft+1):(nfft-1)]中间点(t=0)的相位,而不是2N-1序列[t=0:(2*N-2)]第一点的相位。

补充内容 (2021-11-20 22:01):
楼主认为图2显示的相位46度不正确。图2是t=0:2N-2吋的相位图,它显示的不是初相100度,而是序列第N点的相位46度。可以计算出当初相100时的第N点相46度。

补充内容 (2021-11-20 22:09):
楼主认为图2的相位图显示的46度不正确。图2是t=0:2*N-2时的apfft相位图, 它显示这个序列第N点的相位值是46度。这个序列的第一个样点是100度,笫N点样点的相位就是46度,没有錯。

补充内容 (2021-11-20 22:11):
没有錯

补充内容 (2021-11-21 21:35):
为什么”这个序列的第一个样点相位是100度,第N个相位是46度”
請阅读下面的贴子。有详细的计算过程。
回复此楼

wx_KoU55144 发表于 2021-11-16 20:28:32
你好,我跟你遇到一样的问题了,请问你现在解决了吗?怎么才能改成从0开始模拟实际采样?

wx_KoU55144 发表于 2021-11-16 20:30:12
zhwang554 发表于 2021-7-13 03:42
FFT相位谱中的相位值是序到(t=0:N-1)第一个采样点(t=0)的余弦相位值,即初相位。即x=0.5*cos(2*pi*17.92*t/s ...

王老师您好,请问有办法改成从0开始模拟实际采样吗?能否给个提示

zhwang554 发表于 2021-11-20 01:07:38

第N个样点的相位是46度的详细计算过程。

本帖最后由 zhwang554 于 2021-12-6 04:14 编辑


第2图apfft相位谱中显示的(X=19,Y=46)是指2N-1序列
X=0.5*cos(2*pi*17.92*t/256+100*pi/180)
t=0:2N-2时第N个样点的相位46度。

2N-1序列X的第一个样点的相位是100度,第N个样点的相位可计算如下:
序列中相邻样点的相位增加值都是相同的,相隔一个样点的相位差(度)等於:2*pi*17.92*1/256*360/2/pi=25.2度。
第1个样点相位:100度
第2个:               100+25.2=125.2度
第3个:               100+2*25.2=150.4度

第N个(N=256):mod(100+255*25.2,360)=46度
计算值和图2的测量值一致。

所以,apfft相位谱显示的是2N-1序列中间序列的相位,不管你设置
从0开始的t=0:2N-2, Ljgupup的第2图apfft显示其相位谱,显示当第一个相位是初相100度时,中间第N个样点相位46度。
或 t=-N+1:N-1, Ljgupup的第4图apfft显示其相位谱,显示原点在中间时的初相位100度。
二者都是正确的。但是第4图更好,它将初相100度完美显示,第2图显示的中间样点相位46度,还要计算一下才知对不对。
当然,也可以从中间样点相位值46度,反算出(46度减去255个25.2度)得到的第一个样点的初相位100度是否正确。

当2N-1序列的第一个样点的相位是100度,第N个样点的相位可直接计算如下:
y=mod(2*3.1416*17.92*255/256*360/3.1415/2+100,360)=46度。
计算值和图二的测量值一致,

apfft相位谱显示的是2N-1序列中间样点的相位,不管你设置
t=0:2N-2, t=0,100度, Ly gupup的第2个图apfft显示相位谱而
t=-N+1:N-1, t=0,100度, Lgupup的第4个图apfft显示相位谱,
二者都是正确的,但第4图更㚥它将初相100度完t美显示。第二图显示的l相位还要计算一下才知对不对。
当然,也可以中间样点相位值46度,反算出初相100度。













补充内容 (2022-1-4 01:42):
用matlab计算很简单
先求第N样点的余弦值, 再求它的反余弦度数值
a=x(256)=0.6947,
acosd(a)=46度

补充内容 (2022-1-4 05:58):
用mod(2*pi*17.92*255/256*360/2*pi+100,360)=46度,这个公式求任意样点相位有一个问题,就是当计算得到的相位值在0到180度吋可直接使用。但是当计算的相位在180度到360度时,还要換算一下。详细分析請阅读下面的贴子。

zhwang554 发表于 2021-11-21 22:50:31

从零开始的设置apfft相位谱测到的是笫N个点的相位

本帖最后由 zhwang554 于 2021-12-6 05:37 编辑
wx_KoU55144 发表于 2021-11-16 20:30
王老师您好,请问有办法改成从0开始模拟实际采样吗?能否给个提示

.
Apfft在实用时,只能用从零开始设置,t=0:2*N-1,从相位谱上测到的是序列中间点的相位。到底对不对一时看不出来。
而同样的序列,用FFT测试得到的是第一点的相位,这点FFT是否更好呢?
实际中,从現㘯采集到的一段数据序列,其第一个样点的相位也是不知道的,到底对不对一时也看不出来。
但是我们通过计算机产生的序列对FFT和apfft模拟测试,知道FFT的相位谙测得是初相,apfft测得是第N个点的相位。而且它的精度如何,噪声性能如何,计算机模拟中都搞得很清楚。
所以在用apfft分析信号时,我们测到虽然是第N点相位,加上测到的频率和振幅,这个信号全清楚了,FFT做到的apfft都能做到,在有些性能还比FFT好些。
所以我们不必为apfft测试用t=-N+1:N-1有疑问。好㑰时间为负的采样得不到,apfft不实用,没法用。这个问题经常有人提出,十分简单的就將apfft否决了。
实际上用负时间序列只是测试用,验证apfft的性能,如非周期采样时,apfft测得的相位不变,apfft的水平相位特性等。
实用时,用从零开始的设置,你只要知道相位谱测到的是笫N个点的相位,一切按照常規分析信号。
apfft提出后,有数百篇实用apfft论文,其中许多是用实际采集的数据,它们使用的都是从零开始的设置。


补充内容 (2022-1-1 05:52):
䃼充:好像时间为负采样得不到,apfft不实用没法用。这个问题经常有人提出,十分简单就将apfft否定了。
这网站就有一个.
参見下贴子链接:
https://www.ilovematlab.cn/threa ... _tid=RelatedContent
发表于 2015-4-28 19:46:27
因为这个方法一开始就不对 只能对于对称的构造信号 t=0是关键

补充内容 (2022-1-1 06:18):
楼主比较客观。他说,“哎,真想不明白呀,图2正确说明算法正确啊”。也对箅法正确性产生怀疑。但见到图4又是正确的,又对自己的怀疑抱不定主意,将怀疑提了出来。图2的46度到底是什么?这次这个问题介决很好。apfft提出后,还是首次详细计算出当初相是100度时,第N个相位是46度,

补充内容 (2022-1-22 20:58):
这个方法只适用於对称的信号,t=0是关键”这个观点很正确,还很深刻。等距取样的三角函数的相位值是一条直线,对任意样点都是对称的。全相fft相位谱有平坦特性就是因取样信号相位的对称性质。所以apfft适用於等距取样旳三角函数。不适合非等距取样信号。而仼意函数可分解於三角函数之和,所以apfft适合分析等距取样旳任意函数。

zhwang554 发表于 2021-11-21 22:54:12

从0开始模拟实际采样实例:

本帖最后由 zhwang554 于 2021-12-6 03:13 编辑

wx_KoU55144 发表于 2021-11-16 20:30
王老师您好,请问有办法改成从0开始模拟实际采样吗?能否给个提示

从0开始模拟实际采样实例:

滚动轴承外圈故障信号数据 data130.mat http://home.vibunion.com/blog-62061-17809.html
悬臂梁的脉冲响应数据 采样数据.txt  http://home.vibunion.com/blog-62061-17754.html
电力系统故障信号 i11.mat  http://home.vibunion.com/blog-62061-17752.html
消除肌电信号中的工频信号  http://home.vibunion.com/blog-62061-17746.html
脉动压力数据 x.dat 的功率谱   http://home.vibunion.com/blog-62061-17745.html
求这组数据的初相位??   http://home.vibunion.com/blog-62061-18059.html
如何求幅值  http://home.vibunion.com/blog-62061-18057.html
数据去噪   http://home.vibunion.com/blog-62061-18055.html
滚动体点蚀数据  http://home.vibunion.com/blog-62061-18048.html
测3相交流电数据data1.txt的相位    http://home.vibunion.com/blog-62061-17736.html
Blind Timing Skew Estimation Based on Spectra Sparsity and All Phase FFT for Time-Interleaved ADCs   http://home.vibunion.com/blog-62061-18579.html
百篇全相位应用论文精选  http://home.vibunion.com/blog-62061-18290.html
全相位索引  http://home.vibunion.com/blog-62061-18754.html






zhwang554 发表于 2021-12-5 18:20:07

楼主提供的这个序列是周期取样序列

本帖最后由 zhwang554 于 2021-12-6 08:41 编辑

楼主提供的这个   频率 f=17.92, 取样频率sf=256, 序列长度nfft=256, 初相=100    序列。
仔细分析这个序列每一个样点的相位:
第  56个样点,,
第156个样点,
第256个样点
的相位都是46度。它是一个周期100的周期取样序列。也就是
第1个样点,第101个样点,201个样点的柦位都是100度。它是一个周期序列。
判断一个序列是否周期序列的条件是
f*nfft/sf  是整数。
今  f=17.92,sf=256, nfft=100,
f*N/sf=17.92*100/256=7。
所以这个序列用DFT分析。当序列长度nfft选择为100时,是周期取样序列,用DFT分析得到的相位谱显示正确的初相100度,不须校正。但当nfft不选为100,选择256时,FFT相位谱不显示100度,须测出频率再校正。麻烦多了。
但是选nfft=100,不是2的次方, 不能用FFT快速运算。计算量大幅增加。
当用apfft分析时,不管nfft取什么正整数。用从0开始的设置,相位谱显示46度;用0在中间的设置,相位谱显示100度。注意,不管nfft取任意正整数。所以apfft分析时,可选nfft=256非正周期取样,能够用FFT快速运算。这是apfft的优点: 非正周期取样时,相位值不变,不须校正。所谓的相位不变性。
在实用FFT的仪器中,都须保证周期采样,都须用晶体稳定频率以保证周期采样。但是有的测量不可能做到周期采样,如电力谐波,间谐波测试,谐波还可做到周期采样,间谐波不可能。

Ljgupup 发表于 2021-12-7 09:29:27
zhwang554 发表于 2021-12-5 18:20
楼主提供的这个   频率 f=17.92, 取样频率sf=256, 序列长度nfft=256, 初相=100    序列。
仔细分析这个序列 ...

你好,帮我看看我的另外一个帖子,也是关于全相位的,我用老师的相位校正算法时发现,采样率会对校正精度产生很大影响,不知道我的测试结果正不正确
帖子链接如下:
https://www.ilovematlab.cn/thread-612232-1-1.html

zhwang554 发表于 2022-1-4 05:02:28
zhwang554 发表于 2021-11-20 01:07

第2图apfft相位谱中显示的(X=19,Y=46)是指2N-1序列
X=0.5*cos(2*pi*17.92*t/256+100*pi/180)

求信号x=cos(2*pl*17.92*t/sf+100*pi/180)第256个样点相位

用mod(2*pi*17.92*255/256*360/2*pi+100,360)=46度,这个公式求任意样点相位有一个问题,就是当计算得到的相位值在0到180度吋可直接使用。但是当计算的相位在180度到360度时,还要換算一下。
这是因为余弦信号x的256个取样点的波形如下图。
余弦信号的波形图在0度至180度之间。大於180度还要換祘一下。我们详细分析这个余弦信号x这个实例。
第1个样点 100度。     下一个样点加25.2度,直至180庋
第2个样点 100+25·2=125.2度
第3个样点125.2+25.2=150.4度
第4个样点150.4+25.2=175.6 度
第5个样点175·6+25.2=200.8度。它超过180度20.8度,余弦信号180度以后垂直度数変小,180-20.8=159.2度,图中显示159.2度
第6个样点159.2 - 25.2=134度。   从这点开始每样点相位减25.2直至0度
第7个样点134-25.2=108度
.........
    反复至256个样点
从第5个样点的換算过程可見換算过程还是挺麻烦的。

如果我们用先求a=x(256)=0.6947,
再acosd(a)=46度
这个先求余弦cos值,再求反余弦acosd度数,将取mod和上面的換算步骤都省略了,直接将180度到360度的计算值換算到小於180度。
利用这个方法可以将从第一个样点到第256个样点的相位全部算出如下:

看上去很麻烦
用matlab计算很简单
先算出256个样点的余弦值x,在-1到1之间,再用acosd(x)求反余弦得相位度数,

zhwang554 发表于 2022-1-4 05:33:40
zhwang554 发表于 2022-1-4 05:02
求信号x=cos(2*pl*17.92*t/sf+100*pi/180)第256个样点相位

用mod(2*pi*17.92*255/256*360/2*pi+100,360) ...

信号x=cos(2*pl*17.92*t/sf+100*pi/180)的256个取样点的波形如下图。
b000.jpeg
利用acosd(x)这个方法可以将从第一个样点到第256个样点的相位全部算出如下:
X= 1.0000      2.0000     3.0000      4.0000     5.0000      6.0000     7.0000       8.0000     9.0000    10.0000   11.0000  12.0000   13.0000   14.0000
Y= 100.0000 125.2000  150.4000  175.6000  159.2000  134.0000  108.8000   83.6000   58.4000   33.2000    8.0000   17.2000   42.4000   67.6000

   15.0000   16.0000   17.0000    18.0000    19.0000    20.0000   21.0000    22.0000   23.0000   24.0000   25.0000   26.0000   27.0000   28.0000
   92.8000  118.0000  143.2000  168.4000  166.4000  141.2000  116.0000  90.8000   65.6000   40.4000   15.2000   10.0000   35.2000   60.4000

   29.0000   30.0000   31.0000    32.0000    33.0000   34.0000    35.0000    36.0000   37.0000   38.0000   39.0000   40.0000   41.0000   42.0000
   85.6000  110.8000  136.0000  161.2000  173.6000  148.4000  123.2000  98.0000   72.8000   47.6000   22.4000    2.8000    28.0000   53.2000

   43.0000   44.0000   45.0000    46.0000    47.0000    48.0000    49.0000     50.0000   51.0000   52.0000   53.0000   54.0000   55.0000   56.0000
   78.4000  103.6000  128.8000  154.0000  179.2000  155.6000  130.4000  105.2000   80.0000   54.8000   29.6000     4.4000    20.8000   46.0000

   57.0000   58.0000   59.0000   60.0000    61.0000   62.0000     63.0000    64.0000   65.0000   66.0000   67.0000   68.0000   69.0000   70.0000
   71.2000   96.4000  121.6000  146.8000  172.0000  162.8000  137.6000  112.4000   87.2000   62.0000   36.8000   11.6000   13.6000   38.8000

   71.0000   72.0000   73.0000   74.0000   75.0000     76.0000    77.0000     78.0000    79.0000   80.0000   81.0000   82.0000   83.0000   84.0000
   64.0000   89.2000  114.4000  139.6000  164.8000  170.0000  144.8000   119.6000   94.4000   69.2000   44.0000   18.8000    6.4000   31.6000

   85.0000   86.0000   87.0000   88.0000    89.0000   90.0000    91.0000    92.0000    93.0000   94.0000   95.0000   96.0000   97.0000   98.0000
   56.8000   82.0000  107.2000  132.4000  157.6000  177.2000  152.0000  126.8000  101.6000   76.4000   51.2000   26.0000    0.8000   24.4000

   99.0000  100.0000  101.0000  102.0000  103.0000  104.0000  105.0000   106.0000  107.0000  108.0000  109.0000  110.0000  111.0000  112.0000
   49.6000   74.8000   100.0000  125.2000  150.4000  175.6000  159.2000   134.0000  108.8000   83.6000   58.4000   33.2000     8.0000    17.2000

  113.0000  114.0000  115.0000  116.0000  117.0000  118.0000  119.0000  120.0000  121.0000  122.0000  123.0000  124.0000  125.0000  126.0000
    42.4000   67.6000   92.8000   118.0000  143.2000  168.4000  166.4000   141.2000  116.0000   90.8000   65.6000   40.4000    15.2000    10.0000

  127.0000  128.0000  129.0000  130.0000  131.0000  132.0000  133.0000  134.0000  135.0000  136.0000  137.0000  138.0000  139.0000  140.0000
   35.2000    60.4000   85.6000   110.8000  136.0000  161.2000  173.6000  148.4000  123.2000   98.0000   72.8000    47.6000    22.4000    2.8000

  141.0000  142.0000  143.0000  144.0000  145.0000  146.0000  147.0000  148.0000  149.0000  150.0000  151.0000  152.0000  153.0000  154.0000
   28.0000   53.2000     78.4000  103.6000  128.8000  154.0000  179.2000  155.6000  130.4000  105.2000   80.0000    54.8000   29.6000    4.4000

  155.0000  156.0000 157.0000  158.0000  159.0000  160.0000  161.0000  162.0000  163.0000  164.0000  165.0000  166.0000  167.0000  168.0000
   20.8000    46.0000   71.2000    96.4000   121.6000  146.8000  172.0000  162.8000  137.6000  112.4000    87.2000   62.0000   36.8000    11.6000

  169.0000  170.0000  171.0000  172.0000  173.0000  174.0000  175.0000  176.0000  177.0000  178.0000  179.0000  180.0000  181.0000  182.0000
    13.6000   38.8000   64.0000    89.2000   114.4000  139.6000  164.8000  170.0000  144.8000  119.6000   94.4000    69.2000   44.0000   18.8000

  183.0000  184.0000  185.0000  186.0000  187.0000  188.0000  189.0000  190.0000  191.0000  192.0000  193.0000  194.0000  195.0000  196.0000
     6.4000   31.6000    56.8000   82.0000   107.2000  132.4000  157.6000   177.2000  152.0000  126.8000  101.6000   76.4000    51.2000    26.0000

  197.0000  198.0000  199.0000  200.0000  201.0000 202.0000  203.0000  204.0000  205.0000  206.0000  207.0000  208.0000  209.0000  210.0000
    0.8000     24.4000   49.6000   74.8000   100.0000  125.2000  150.4000  175.6000  159.2000  134.0000  108.8000   83.6000    58.4000   33.2000

  211.0000  212.0000  213.0000  214.0000  215.0000  216.0000  217.0000  218.0000  219.0000  220.0000  221.0000  222.0000  223.0000  224.0000
     8.0000   17.2000    42.4000   67.6000    92.8000    118.0000  143.2000 168.4000  166.4000  141.2000   116.0000    90.8000   65.6000   40.4000

  225.0000  226.0000  227.0000  228.0000  229.0000  230.0000  231.0000  232.0000  233.0000  234.0000  235.0000  236.0000  237.0000  238.0000
    15.2000   10.0000   35.2000    60.4000    85.6000  110.8000  136.0000  161.2000  173.6000  148.4000  123.2000    98.0000   72.8000   47.6000

  239.0000  240.0000  241.0000  242.0000  243.0000  244.0000  245.0000  246.0000  247.0000  248.0000  249.0000  250.0000  251.0000  252.0000
    22.4000    2.8000     28.0000   53.2000   78.4000   103.6000  128.8000  154.0000  179.2000  155.6000  130.4000  105.2000   80.0000   54.8000

  253.0000  254.0000  255.0000  256.0000
    29.6000    4.4000   20.8000   46.0000


zhwang554 发表于 2022-1-7 02:43:43
zhwang554 发表于 2022-1-4 05:33
信号x=cos(2*pl*17.92*t/sf+100*pi/180)的256个取样点的波形如下图。

利用acosd(x)这个方法可以将从第一 ...


上贴子所示信号相位图虽显示了256个样点的相位,但是它的连线在靠近0度和180度附近图形畸変。下面的程序显示改进后的正确图形。
b000good111new.jpeg
clear
clc
close all hidden
%准备数据
sf=256*1; %采样频率
nfft=256*1; %FFT计算点数
n=2*nfft-1; %数据采样点数

%采样点
t=(-nfft+1):(nfft-1);
t=0:1:(2*nfft-2);
tt=1:(2*nfft-1);

x=1*cos(2*pi*17.92*t/sf+100*pi/180); %采样函数频率17.92,相位100°

%相位值为180度的t值
a=(180-100+360*(0:17))/25.2+1;
%相位值为0度的t值
b=(-100+360*(1:18))/25.2+1;
c=[1   a(1) b(1) a(2) b(2) a(3) b(3) a(4) b(4) a(5) b(5) a(6) b(6) a(7) b(7) a(8) b(8) a(9) b(9) a(10) b(10) a(11) b(11) a(12) b(12) a(13) b(13) a(14) b(14) a(15) b(15) a(16) b(16) a(17) b(17) a(18) b(18) 256 ];
d=[100 180  0    180   0   180  0     180  0   180   0   180   0    180  0    180   0   180  0     180  0     180   0     180   0     180    0    180   0     180    0     180    0   180   0   180   0  46];
plot(c,d,'x')
hold on
plot(c,d,'r')
hold on
y=acosd(x);
plot(y(1:256),'x')
title('cos(2*pi*17.92*(1:256)/256+100*pi/180 phase diagram');
xlim([0 256]);ylim([0 180]);
xlabel('samples');ylabel('degree');
hold on


cosphase.m (1.03 KB, 下载次数: 1)

zhwang554 发表于 2022-1-7 07:23:00
zhwang554 发表于 2022-1-7 02:43
上贴子所示信号相位图虽显示了256个样点的相位,但是它的连线在靠近0度和180度附近图形畸変。下面的程序 ...

上贴画出的余弦取样信号相位图, 由于余弦信号0度在上方, 180在下方, 所以上图不是matlab程序cosphase.m直接画出,而要在paint制图软件上作水平軸反转才能显示0度在上方,180度在下方的余弦取样相位图。由于图形不是matlab直接画出,虽然图中“x”表示取得点,但是不能在图中显示取样点x为多少,相位y为多少。必须对照acodd(x)相位读数才知相位值。

下面给出正弦取样相位图,相位值在+90度和 -90之间,0軸在中间。在matlab中可直接画出。matlab直接画出的画中每个样点的x和Y相位值可以在画中显示。十分方便,直观。在正弦相位图中可以清楚显示任何一个样点的相位值。 sinphase.jpg
clear
clc
close all hidden
%准备数据
sf=256*1; %采样频率
nfft=256*1; %FFT计算点数
n=2*nfft-1; %数据采样点数

%采样点
t=(-nfft+1):(nfft-1);
t=0:1:(2*nfft-2);
tt=1:(2*nfft-1);

x=1*sin(2*pi*17.92*t/sf+0*pi/180); %采样函数频率17.92,相位100°

a=(90-0+360*(0:18))/25.2+1;
b=(-90-0+360*(1:18))/25.2+1;
c=[1  a(1) b(1) a(2) b(2) a(3) b(3) a(4) b(4) a(5) b(5) a(6) b(6) a(7) b(7) a(8) b(8) a(9) b(9) a(10) b(10) a(11) b(11) a(12) b(12) a(13) b(13) a(14) b(14) a(15) b(15) a(16) b(16) a(17) b(17) a(18) b(18) a(19) 256 ];
d=[0 90  -90    90  -90   90  -90   90  -90   90  -90   90  -90  90   -90   90   -90  90   -90   90   -90   90   -90   90    -90    90   -90    90    -90    90   -90    90   -90   90   -90    90   -90  +90    -54];
plot(c,d,'x')
hold on
plot(c,d,'r')
hold on
y=asind(x);
plot(y(1:256),'x')
title('sin(2*pi*17.92*(1:256)/256+0*pi/180 phase diagram');
xlim([1 256]);ylim([-90 90]);
xlabel('samples');ylabel('degree');
hold on




zhwang554 发表于 2022-1-7 23:26:07
zhwang554 发表于 2022-1-7 07:23
上贴画出的余弦取样信号相位图, 由于余弦信号0度在上方, 180在下方, 所以上图不是matlab程序cosphase.m直 ...

上贴画出的正弦取样的相位图,它呈现直线型。
本贴再配上正弦取样的图形,它呈正弦曲线型。
两者配合对照看,可见每一样点在正弦曲线上的位置,正弦值和相位值。 sinphasesamplesnew.jpeg
sf=256; %采样频率
nfft=256; %FFT计算点数
n=2*nfft-1; %数据采样点数

%采样点
t=(-nfft+1):(nfft-1);
t=0:1:(2*nfft-2);
tt=1:(2*nfft-1);

x=1*sin(2*pi*17.92*t/sf+0*pi/180); %采样函数频率17.92,相位100°

a=(90-0+360*(0:18))/25.2+1;
b=(-90-0+360*(1:18))/25.2+1;
c=[1  a(1) b(1) a(2) b(2) a(3) b(3) a(4) b(4) a(5) b(5) a(6) b(6) a(7) b(7) a(8) b(8) a(9) b(9) a(10) b(10) a(11) b(11) a(12) b(12) a(13) b(13) a(14) b(14) a(15) b(15) a(16) b(16) a(17) b(17) a(18) b(18) a(19) 256 ];
d=[0 90  -90    90  -90   90  -90   90  -90   90  -90   90  -90  90   -90   90   -90  90   -90   90   -90   90   -90   90    -90    90   -90    90    -90    90   -90    90   -90   90   -90    90   -90  +90    -54];
subplot(2,1,2)
plot(c,d)
hold on
plot(c,d,'r')
hold on
y=asind(x);
plot(y(1:256),'*')
title('sin(2*pi*17.92*(1:256)/256+0*pi/180) samples phase diagram');
xlim([1 256]);ylim([-90 90]);
xlabel('samples');ylabel('degree');
hold on
subplot(211)
plot(x,'x')
title('sin(2*pi*17.92*(1:256)/256+0*pi/180) samples diagram');
xlim([1,256])
hold on
plot(x,'r')
xlabel('samples');ylabel('degree');


sinphase1.m (1.15 KB, 下载次数: 1)

zhwang554 发表于 2022-4-22 01:58:27
当一个频率成份振幅是很小的时,很易被泄漏淹没, 这个频率成分在FFT和apfft振幅谱不显示。但这个频率成分在apfft相位谱显示的是采样第一个样点的余弦初相,不管这个频率成分的振幅多小,其初相是0到360度之间的一个数傇,取中间偏小值可能也有几十度,其显示效果远远超过很小的振幅值。
apfft相位谱具有水平相位特性,即在振幅峰值附近一定频率范围內相位值相同,峯值附近的相位值是峰值矢量泄漏出去的,其实数和虚数部份都很小,所以可以加一个阀值,使峰值显示相位值,周同围相位值置零。加阈值的方法见附程序个
这样apfft的相位谱更明显显示某一成分的相位值
通过apfft所得出的结果的精度十分可观。
异常眀显的相位変化使得该方法的效果远远超过振幅谱,
%全相傅里叶变换
clear;clc;close all hidden
%准备数据
sf=256; %采样频率
nfft=256; %FFT计算点数
n=2*nfft-1; %数据采样点数
t=(-nfft+1):(nfft-1);
x=cos(50.3*2*pi*t/sf+50*pi/180)+0.001*cos(85.3*2*pi*t/sf+100*pi/180);%+cos(121.4*2*pi*t/sf+pi/2);
win1=hanning(nfft)';
win2=hanning(nfft)';
win=conv(win1,win2); %卷积窗
win=win/sum(win);     %归一化
w=pi*2;
y1=x.*win;                %加权
y1a=y1(nfft:end) + [0 y1(1:nfft-1)];
Y2=fft(y1a, nfft); %apFFT
s=0.0001;
q=real(Y2);
for i=1:256;
if q(i)<s;
    q(i)=0;
elseif  q(i)>s;
    q(i)=q(i);
end
end
qq=imag(Y2);
for i=1:256;
if qq(i)<s;
    qq(i)=0;
elseif  qq(i)>s;
    qq(i)=qq(i);
end
end
Y2=q+j*qq;
Y1=fft(x(1:nfft).*win1,nfft); %传统FFT
s=0.1;

a1=abs(Y1);
a2=abs(Y2);
p1=angle(Y1)*180/pi;
p2=angle(Y2)*180/pi;

subplot(2,2,1);
plot(a1)
xlim([0 256/2]);
title('传统FFT');
subplot(2,2,2);
plot((a2)); %幅度谱
xlim([0 256/2]);
title('apFFT');
subplot(2,2,3);
plot(p1-s); %相位谱
xlim([0 256/2]);
subplot(2,2,4);
plot(p2); %相位谱
xlim([0 256/2]);

zhwang554 发表于 2022-4-22 02:16:27
zhwang554 发表于 2022-4-22 01:58
当一个频率成份振幅是很小的时,很易被泄漏淹没, 这个频率成分在FFT和apfft振幅谱不显示。但这个频率成分在 ...

aaaaa.jpg

zhwang554 发表于 2022-5-1 20:45:36
zhwang554 发表于 2022-4-22 01:58
当一个频率成份振幅是很小的时,很易被泄漏淹没, 这个频率成分在FFT和apfft振幅谱不显示。但这个频率成分在 ...

FFT和apfft谱分析比较程序
常见的fft和apfft比较程序中,apfft相位谱呈现水平相位特性。在上贴中,将apfft相位谱设限,相位谱只显示振幅谱峰值频率成分的相位,这样当这个频率成分在掁幅值很小振幅谱显示不清楚时,相位谱显示的初相值仍较大,能够清晰显示这个频率成分的存在。但是用阈值的程序不是很好,因为须增加阈值的选择。本贴中采用捜索apfft振幅谱的全部峰值的频率,再求各峰值频率的相位值。这样在apfft的相位谱中清楚正确显示各成分的初相。同样方法也用於fft相位谱,先求出fft掁谱峰值的频率,再求峰值频率处的相位值 ,这样fft的相位谱也不再是乱七八糟的,但是它显示的相位值不是初相,还须校正。clear;clc;close all hidden
Fs=256*2;
nfft=256*2;
n=2*nfft-1;
t=(-nfft+1):(nfft-1);
x=100*cos(2*pi*t*50.2/Fs+50*pi/180)+2*cos(2*pi*t*100.3/Fs+100*pi/180)+3*cos(2*pi*t*150.1/Fs+150*pi/180)+4*cos(2*pi*t*180.3/Fs+200*pi/180)+5*cos(2*pi*t*200.3/Fs+140*pi/180);
win=hanning(nfft)';
win1=conv(win,win);
win2=win1/sum(win1);     
y1=x.*win2;               
y1a=y1(nfft:end) + [0 y1(1:nfft-1)];
F1=fft(y1a, nfft);
a1=abs(F1);
p1=angle(F1);
L=fix(nfft/2)+1;
f=(0:L)*Fs/nfft;
for i=1:5
[A,Ind(i)]=max(a1(1:L));
k=Ind(i);
peak_F1(i)=mod(angle(F1(k))*180/pi,360);
StartP=Ind(i)-2;
StopP=Ind(i)+2;
if Ind(i)-2<1
StartP=1;
end
if Ind(i)+2>L
StopP=L;
end
Indx=StartP:StopP;
a1(Indx)=zeros(1,length(Indx));
end
subplot(2,2,2);
a1=abs(F1);
plot(a1); xlim([0 nfft/2]);ylim([0 60]);
title('apFFT amplitude spectrum');
subplot(2,2,4);
stem(Ind, peak_F1,'*r');xlim([0 nfft/2]);ylim([0 360]);
title('apFFT phase spectrum');

Fs=256*2; %采样频率
nfft=256*2; %FFT计算点数
n=2*nfft-1; %数据采样点数
t=(-nfft+1):(nfft-1);
x=100*cos(2*pi*t*50.2/Fs+50*pi/180)+2*cos(2*pi*t*100.3/Fs+100*pi/180)+3*cos(2*pi*t*150.1/Fs+150*pi/180)+4*cos(2*pi*t*180.3/Fs+200*pi/180)+5*cos(2*pi*t*200.3/Fs+140*pi/180);
win1=hanning(nfft)';
winn=win1/sum(win1);
x1=x(1:nfft);
y1=x1.*winn;                %加权
F2=fft(y1); %传统FFT
a2=abs(F2);
p2=angle(F2);
L=fix(nfft/2)+1;
f=(0:L)*Fs/nfft;
for i=1:5
[A,Ind(i)]=max(a2(1:L));
k=Ind(i);
peak_F2(i)=mod(angle(F2(k))*180/pi,360);
StartP=Ind(i)-2;
StopP=Ind(i)+2;
if Ind(i)-2<1
StartP=1;
end
if Ind(i)+2>L
StopP=L;
end
Indx=StartP:StopP;
a2(Indx)=zeros(1,length(Indx));
end
a2=abs(F2);
subplot(2,2,1);
plot(a2);
title('FFTamplitude spectrm');
xlim([0 nfft/2]);ylim([0 60]);
subplot(2,2,3);
stem(Ind, peak_F2,'*r');
title('FFT phase spectrum');
xlim([0 nfft/2]);ylim([0 360]);

fft_apfft.jpg

zhwang554 发表于 2022-5-6 04:25:41
zhwang554 发表于 2022-5-1 20:45
FFT和apfft谱分析比较程序
常见的fft和apfft比较程序中,apfft相位谱呈现水平相位特性。在上贴中,将apff ...

FFT和apfft谱分析比较程序


上贴的apfft相位谱只显示振幅峰值处频率的相位值,和常用的apfft水平相位谱不一样,这样能清晰显示振幅值较小的频率成分,虽然其在振幅谱中显示不明显,但在apfft相位谱中显示明显。但apfft相位谱还有一个性能,即振幅峰值对应的相位谱如果是斜线,则这个振幅波不是单一频率成分。这个很好的性能在上贴程序中没有反映出来。
本贴程序介决了这个问题即对每一个振幅峰值处的频率值k,显示k,k+1,k-1三处的相位值,若峰值对应的的三个相位值是水平的,这个峰值对应的是单一频率成份若三个相位值是水平的,则不是单一频率成分。
如下程序中有四个频率成分,其中111.9 和111.2是在一个振幅峰值内,所以对应的三个相位谱不是水平的。
在另一附图中显示二种apfft相位谱对照图。
c20220505b.jpg c20220505a.jpg clear;clc;close all hidden
Fs=256*2;
nfft=256*2;
n=2*nfft-1;
t=(-nfft+1):(nfft-1);
x=100*cos(2*pi*t*50.0/Fs+100*pi/180)+49*cos(2*pi*t*111.9/Fs+120*pi/180)+18*cos(2*pi*t*112.2/Fs+40*pi/180)+7*cos(2*pi*t*160.4/Fs+170*pi/180);
win=hanning(nfft)';
win=hanning(nfft)';
win1=conv(win,win);
win2=win1/sum(win1);     
y1=x.*win2;               
y1a=y1(nfft:end) + [0 y1(1:nfft-1)];
F1=fft(y1a, nfft);
a1=abs(F1);
p1=angle(F1);
L=fix(nfft/2)+1;
f=(0:L)*Fs/nfft;
for i=1:3;
win=hanning(nfft)';
[A,Ind(i)]=max(a1(1:L));
k=Ind(i);
peak_F1(i)=mod(angle(F1(k))*180/pi,360);
peak_F11(i)=mod(angle(F1(k+1))*180/pi,360);
peak_F111(i)=mod(angle(F1(k-1))*180/pi,360);
StartP=Ind(i)-2;
StopP=Ind(i)+2;
if Ind(i)-2<1
StartP=1;
end
if Ind(i)+2>L
StopP=L;
end
Indx=StartP:StopP;
a1(Indx)=zeros(1,length(Indx));
end
subplot(2,2,2);
a1=abs(F1);
Ind2=Ind+1;
Ind3=Ind-1;
plot(a1); xlim([0 200]);ylim([0 60]);
title('apFFT amplitude spectrum');
subplot(2,2,4);
stem(Ind, peak_F1,'_r');
xlim([0 200]);ylim([0 360]);
title('apFFT phase spectrum');
hold on
stem(Ind2, peak_F11,'_r');xlim([0 200]);ylim([0 360]);
hold on
stem(Ind3, peak_F111,'_r');xlim([0 200]);ylim([0 360]);


win1=hanning(nfft)';
winn=win1/sum(win1);
x1=x(nfft:2*nfft-1);
y1=x1.*winn;                %加权
F2=fft(y1); %传统FFT
a2=abs(F2);
p2=angle(F2);
L=fix(nfft/2)+1;
f=(0:L)*Fs/nfft;
for i=1:3
[A,Ind(i)]=max(a2(1:L));
k=Ind(i);
peak_F2(i)=mod(angle(F2(k))*180/pi,360);
StartP=Ind(i)-2;
StopP=Ind(i)+2;
if Ind(i)-2<1
StartP=1;
end
if Ind(i)+2>L
StopP=L;
end
Indx=StartP:StopP;
a2(Indx)=zeros(1,length(Indx));
end
a2=abs(F2);
subplot(2,2,1);
plot(a2);
title('FFTamplitude spectrm');
xlim([0 200]);ylim([0 60]);
subplot(2,2,3)
stem(Ind, peak_F2,'*r');
title('FFT phase spectrum');
xlim([0 200]);ylim([0 360]);

figure
subplot(311)
plot(a1)
xlim([0 200])
title('apFFT amplitude spectrum ');
subplot(312)
stem(Ind, peak_F1,'_r');
xlim([0 200]);ylim([0 360]);
title('apFFT phase spectrum A');
hold on
stem(Ind2, peak_F11,'_r');xlim([0 200]);ylim([0 360]);
hold on
stem(Ind3, peak_F111,'_r');xlim([0 200]);ylim([0 360]);
subplot(313)
plot(mod(p1*180/pi,360),'r')
title('apFFT phase spectrum B');
xlim([0 200])


zhwang554 发表于 2022-5-15 06:38:21
zhwang554 发表于 2022-5-6 04:25
FFT和apfft谱分析比较程序

apfft 水平相位特性细节apfft 相位谱呈水平相位特性,即在振幅谱峰值频率的二旁,柤位值是相同的,都等於初相。
我们研究一个简单实例。N=32,f=5.3,  fs=32, 初相45度。数据予处理后作FFT的输出为复数
y2_fft(1:N/2) =
Columns 1 through 4
  2.8142e-06 + 2.8142e-06i  9.9789e-06 + 9.9789e-06i   5.0221e-05+ 5.0221e-05i   0.00050368 + 0.00050368i
Columns 5 through 8
   0.048261 +   0.048261i    0.62451 +    0.62451i    0.35244 +   0.35244i   0.0026577 +  0.0026577i
Columns 9 through 12
  7.9568e-05 + 7.9568e-05i   8.002e-06 +  8.002e-06i   1.2827e-06 + 1.2827e-06i   2.5347e-07 + 2.5347e-07i
Columns 13 through 16
   5.268e-08 +  5.268e-08i   9.6302e-09 + 9.6302e-09i    1.001e-09 + 1.001e-09i   6.0352e-12 +6.0353e-12i
p2(1:N/2)=
Columns]
]
]
   45      45     45     45     45     45     45.001


由於初相为45度,其实数部分数值和虚数相等,频率k=6处数值为 0.62451 + 0.62451i ,数值最大。k=7处为 0.35244 +  0.35244i,k=9
处为 7.9568e-05 + 7.9568e-05i ,实部和虚部数值迅速缩小,但数值仍相等但有些误差,离k=6越远,实虚数仍相等,但数值更小为e-10.


close all;clc;clear all;
N=32;fs=32;f=5.3;
t=-N+1:N-1;w=2*pi;
y=1.*exp(j*(w*t*f/fs+pi/4));
y1 = y(N:2*N-1);
win = hanning(N)';
win1 = win/sum(win);
y11= y1.*win1;
y11_fft = fft(y11,N);
a1 = abs(y11_fft);
p1 = mod(angle(y11_fft)*180/pi,360);
y2 = y(1:2*N-1);
winn = conv(win,win);
win2 = winn/sum(winn);
y22= y2.*win2;
y222=y22(N:end)+[0 y22(1:N-1)];
y2_fft = fft(y222,N);
y2_fft(1:N/2)
a2 = abs(y2_fft);
p2=mod( angle(y2_fft)*180/pi,360);
p2(1:N/2)
tt=1:N;
subplot(211);plot(tt,10*log10(a1),'b',tt,10*log10(a2),'r');
subplot(211);plot(tt,a1,'b',tt,a2,'r');
title('(b) 对数振幅谱  ');
%ylim([1,-12]);
xlim([0,N/2]);
legend('N=32   fft ','N=256   apfft ');
xlabel('f');ylabel('db');
grid on
grid minor
subplot(212);plot(tt,p1,'b',tt,p2,'r');
title('(a) N=32  相位谱 ');
ylim([0,400]);xlim([0,N/2]);
legend('N=32  fft','N=32   apfft ');
xlabel('f');ylabel('度');
grid on
grid minor




您需要登录后才可以回帖 登录 | 注册

本版积分规则

相关帖子
相关文章
热门教程
站长推荐
快速回复 返回顶部 返回列表