如何求幅值？？

kleeblatt

急求问各位，如下情况如何求信号幅值：
已知：信号图形Asin(wt*t+phi) ，基波初相位phi，基波频率w
求：幅值
但算法需在simulink中实现，故不可直接使用Asin(wt+phi)/ sin(wt+phi)，因sin(wt+phi)存在零值。求问各位还能有什么方法实现幅值的输出，谢谢！！！

flyback

这......
直接读最大值不就可以了

kleeblatt

信号Asin(wt+phi)是经dsp处理的，所以读取dsp后的最大值很可能不是原信号的最大值:@o

zhwang554

　　用离散频谱校正可以测出正弦函数的振幅, 有多种校正法都可用,有些校正精度很高. 这些方法大都利用插值公式计算振幅值
　　但有特点的是fft/apfft校正法.它与窗函数无关,和频偏值也无关,由fft和apfft两个振幅值决定.
　　FFT和apFFT组合校正法须2N-1个样点,对前N个样点作FFT,得FFT振幅谱a1,对全部2N-1个样点作apFFT,得apFFT振幅谱a2, FFT功率谱/apFFT振幅谱 , a1.^2/a2 , 即振幅校正谱.
　　程序如下
　　close all;clc;clear all;
　　N=256;
　　t=-N+1:N-1;
　　A=110.1111;f=49.2;phi=40;
　　y=A*cos(2*pi*t*f/N+phi*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个输入数据
　　win =hanning(N)';;
　　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振幅谱
　　aa=(a1.^2)./a2*2;%振幅校正值
　　r=round(f)+1;%振幅谱峰值处频率
　　disp('振幅校正值')
　　aaa=aa(r)
　　运行结果
　　振幅校正值
　　aaa=110.1111
　　[ 本帖最后由 zhwang554 于 2009-11-15 09:18 编辑 ]

kleeblatt

请问a1的值是提前给定的么？我现在只有sin(wt+phi)的值和Asin(wt=phi)的值，如何给定A呢？

zhwang554

抱歉, 原来程序中用了2个a1, 让人弄不清楚,  改了一下

笫一个a1改为A, 即信号Asin(wt=phi)的振幅

第二个a1(保留在程序中)是序列y11的FFT的振幅谱

          其中y11是加hanning窗y1的序列

          其中y1是给定的余弦信号Asin(wt=phi)的取样值

你说: "我现在只有sin(wt+phi)的值和Asin(wt=phi)的值，如何给定A"?
离散频谱校正法只须Asin(wt=phi)的值，就能定A
离散频谱校正是从信号Asin(wt+phi)取样值的DFT的振幅谱测A, 而且是在非整周期取样下测A. 整周期取样下Asin(wt+phi)取样值的DFT的振幅值就是A.

这个方法是离散频谱校正法, 对余弦信号Asin(wt+phi)的取样值作FFT得振幅谱a1, 作apFFT得振幅谱a2,   a1.^2/a2 即得A (在振幅谱a1峰值频率处)

关于离散频谱校正法论坛中贴子很多, 可参看相关贴子,如
http://forum.vibunion.com/UChome/space-62061-do-blog-id-17741.html

[ 本帖最后由 zhwang554 于 2009-11-15 20:51 编辑 ]

kleeblatt

sorry that i can't type chinese at the moment.
want to ask again, that i can't give the function of the input signal, just 128 discrete values in 1 sec. how can i get the apfft values, thanks!!!!!

zhwang554

　　你的问题是:若1秒内得128个采样,如何作apFFT?
　　选定fft阶数N后,当采样频率fs=N赫/秒,作apFFT需2秒采样, 当采样频率fs=N/2赫/秒,作apFFT只需1秒采样,
　　所以若得1秒内128个采样(即采样频率为128赫/秒),相当采样频率为64赫/秒2秒内得128个采样. 计算频率时注意采样频率为64赫/秒.
　　若只得128个采样,你用127采样作N=64的apFFT即可,反正作N阶FFT/apFFT须2*N-1个采样.如何作apFFT 5楼程序中都写出了
　　类同的有一种FFT时移相位差法,它须二组相邻序列作FFT才能校正振幅等,所以也要采样2秒
　　由信号频谱求振幅的一个长处是一次FFT可得多种频率式成份的振幅值. 用fft/apfft校正法得的振幅校正谱可用下aaaa.doc附图示
　　N=256;
　　t=-N+1:N-1;
　　A=1.0;f=29.2;phi=40;
　　y=A*cos(2*pi*t*f/N+phi*pi/180)+2*A*cos(2*pi*t*2*f/N+phi*pi/180)+3*A*cos(2*pi*t*3*f/N+phi*pi/180);
　　附图是3个成分信号y的振幅谱校正图;(a)为fft振幅谱a1,(b)为apfft振幅谱a2,(c)校正振幅谱.
　　3个成分信号振幅校正为1,2,3.而且在峰值附近形成一个平台,很醒目. 有噪时峰值处最准.
　　振幅校正值(111db)
　　0.9999 1.9999 3.0000
　　振幅校正值(10db)
　　1.0536 1.9113 3.0705
　　从图也可见,未经校正的(a)fft振幅谱a1和(b)apfft振幅谱a2峰值都不是1,2,3. 只有校正后(C)的峰值才是1,2,3.
　　 aaaa.doc (30.5 KB, 下载次数: 5)

2009-11-17 21:45 上传

点击文件名下载附件

　　[ 本帖最后由 zhwang554 于 2009-11-17 23:33 编辑 ]

kleeblatt

here is my program:
N=64;
y=phase1;
y1 = y(N:2*N-1);
win =  hanning(N)';;
win1 = win/sum(win);
y11= y1.*win1;
y11_fft = fft(y11,N);
y1_amplitude = abs(y11_fft);
y2 = y(1:2*N-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)];
y2_fft = fft(y222,N);
a2 = abs(y2_fft);
aa=(a1.^2)./a2*2;
r=round(f)+1;
aaa=aa(r);

it shows:

??? Error using ==> times
Matrix dimensions must agree.

Error in ==> Untitled3 at 7
y11= y1.*win1;

kleeblatt

phase1 is an array with 128 sample values in 1 sec.

zhwang554

除红色外,程序对


N=64;
t=-N+1:N;
A=1.0;f=17.2;phi=40;
y=A*cos(2*pi*t*f/N+phi*pi/180);%产生一个128采样数据,可换你的phase1
y1 = y(N:2*N-1);
win =  hanning(N)';;
win1 = win/sum(win);
y11= y1.*win1;
y11_fft = fft(y11,N);
y1_amplitude = abs(y11_fft);
y2 = y(1:2*N-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)];
y2_fft = fft(y222,N);
a2 = abs(y2_fft);
aa=(y1_amplitude.^2)./a2*2;
r=round(f)+1;
aaa=aa(r)

[ 本帖最后由 zhwang554 于 2009-11-18 01:17 编辑 ]

zhwang554

喔!我明白了,不能出现解析式, 只有1秒128个采样(that i can't give the function of the input signal, just 128 discrete values in 1 sec.)

N=64;
phase1=[0.7660    0.0217   -0.7374   -0.9940   -0.5732    0.2382    0.8873    0.9317    0.3412   -0.4818   -0.9765   -0.8058   -0.0859    0.6925    0.9990...
    0.6247   -0.1753   -0.8558   -0.9532   -0.4010    0.4245    0.9607    0.8422    0.1498   -0.6447   -0.9998   -0.6736    0.1116    0.8208    0.9706...
    0.4590   -0.3654   -0.9408   -0.8751   -0.2130    0.5942    0.9965    0.7197   -0.0475   -0.7824   -0.9841   -0.5152    0.3048    0.9171    0.9044...
    0.2754   -0.5413   -0.9891   -0.7629   -0.0168    0.7407    0.9935    0.5692   -0.2430   -0.8896   -0.9299   -0.3366    0.4861    0.9776    0.8028...
    0.0810   -0.6960   -0.9987   -0.6208    0.1801    0.8584    0.9516    0.3964   -0.4290   -0.9620   -0.8395   -0.1449    0.6485    0.9999    0.6699...
   -0.1165   -0.8236   -0.9694   -0.4546    0.3700    0.9425    0.8727    0.2082   -0.5982   -0.9969   -0.7163    0.0525    0.7855    0.9832    0.5109...
   -0.3095   -0.9190   -0.9023   -0.2706    0.5454    0.9898    0.7596    0.0118   -0.7441   -0.9929   -0.5651    0.2478    0.8918    0.9281    0.3319...
   -0.4904   -0.9786   -0.7999   -0.0761    0.6996    0.9985    0.6170   -0.1850   -0.8609   -0.9501   -0.3919    0.4334    0.9634    0.8368    0.1400...
   -0.6522   -1.0000   -0.6663    0.1215    0.8264    0.9682    0.4502   -0.3746];
y=phase1;
y1 = y(N:2*N-1);
win =  hanning(N)';;
win1 = win/sum(win);
y11= y1.*win1;
y11_fft = fft(y11,N);
y1_amplitude = abs(y11_fft);
y2 = y(1:2*N-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)];
y2_fft = fft(y222,N);
a2 = abs(y2_fft);
aa=(y1_amplitude.^2)./a2*2;
[A,r]=max(a2(1:N/2));
aaa=aa(r)

运行结果:
aaa = 1.000

[ 本帖最后由 zhwang554 于 2009-11-18 22:33 编辑 ]

kleeblatt

exactly what i mean!
but i can't see any difference between my program and what u made. pls point ...:@L

kleeblatt

it still shows the same problem:???
Error using ==> times
Matrix dimensions must agree.
phase1 is the signal.
let me describe a little bit about phase1.Actually phase1 should be phase1.signals.values, from 'to workspace'.To workspace's variable name is phase1, and saved format is
structure with time.Sample time is 1/127,so that i get 128 discrete values.
anything wrong with this part?

zhwang554

　　13楼程序在matlab上运行没错.
　　phase1的128个数据由下产生
　　t=0:1/127:1;
　　A=1.0;f=17.2;phi=40;
　　phase1=A*cos(2*pi*t*f+phi*pi/180)
　　你问:"but i can't see any difference between my program and what u made. pls point ..."
　　no difference between your program(10楼) and what i made(13楼).
　　其实, 5楼,10楼,12楼,13楼都是一个程序
　　[ 本帖最后由 zhwang554 于 2009-11-18 22:52 编辑 ]