求助幅值谱图

秋月

请大家帮忙看一下达芬方程的幅值谱图为什么会这样??
那有问题?刚开始还不太懂,请多多指教!
谢谢!
clc
clear all
x0=[0;1.0; 1];
[t,y]=ode45('duffing',[0:0.01:100],x0);
%频谱图
fs=200;          %采样频率
N=10000;           %采样点数
x=y(:,1);  %生成信号（Y轴方向的位移）
Y=fft(x,N);
Pyy = Y.* conj(Y) /N;
df=fs/(N-1);          %频率分辨率Hz
f=(0:N-1)*df;    %频域序列
figure
plot(f,Y);
xlabel('Frequency(HZ)');
ylabel('Magnitude');
title('频率幅值图');
hold on

子程序
function dy=duffing(t,x)
m=1.0;
c=0.1;
k=1.0;
F=13;
w=1.0;
dy=[x(2);
   -c*x(2)-x(1)^3+F*cos(x(3));
   w];

shenyongjun

两个问题：
1、采样频率为100，不是200；
2、你画的是双边谱，一般情况下画单边谱即可

figure
plot(f(1:fix(N/2)),Y(1:fix(N/2)));

秋月

谢谢!不过结果只是少了右边的一半啊,
请高手帮忙看看,指点一下!
谢谢了!

sssssxxxxx921

请注意单边谱和双边谱的表述   这本身就是一个单边谱
另外，你给出的只是采样长度和显示频率的关系，看看采样定理
第三，秋月是绘制的一个功率谱    而非幅值谱

咕噜噜

^_^，申老师可是很厉害的老师啊，虽然这个我不懂，不过申老师的意见还是需要好好参考的

octopussheng

我用我的方法做了一下，你做个参考吧

function df=dafen(t,x,flag,force)
df=[x(2);force*cos(1.2*t)-x(1)^3+x(1)-0.3*x(2)];

clear
ff=0.222;
options=odeset('RelTol',1e-7);
t0=0;
tf=500;
[t,x]=ode45(@dafen,[t0 tf],[0,0],options,[],ff);
plot(x(:,1),x(:,2))

>> Y=fft(xx);                            ％快速fft变换
>> N=length(Y);                       ％求Y的向量长度
>> Y(1)=[];                              ％去掉Y的第一个数据，它是所有数据之和
>> power=abs(Y(1:N/2)).^2;      ％求功率谱
>> nyquist=1/2;
>> freq=(1:N/2)/(N/2)*nyquist;   ％求频率
>> plot(freq,power);grid on;

结果如图示

[ 本帖最后由 octopussheng 于 2007-9-21 09:02 编辑 ]

sssssxxxxx921

小咕  我这可得批评你几句
我们都知道研究生和博士 还是哪个专家他只是在自己的专业领域很精通   并不是说所有的都在行
当然并不是说没有这样的人   也不是说贬低申永君  说实话  我也看过他的论文  关于多频激励的  
但是他是用近似理论求解的  如果我没说错他的老师是杨绍普  石家庄铁道学院  我一开始从图书馆借到的非线性书就有他的一本
在频谱分析方面我不了解他    但我了解我自己   而且我所说的都是非常确信的   例如单边谱和双边谱   功率谱和幅值谱的区别 这些是我非常肯定的，也是我在做程序和看书后获得的
我觉得在上面那几个问题上只有对和不对之分    没有商榷的任何余地  
同时我希望所有的人都尊重权威    不轻易否定权威，但是绝对不要迷信权威
最后声明：我只是就是就事论事   不针对任何人
谢谢  
:victory:

sssssxxxxx921

freq=(1:N/2)/(N/2)*nyquist;   ％求频率
这个采样频率为1时是正确的，但当不是1的时候要注意这个采样频率1/dt  
而且这是还是一个功率谱   和幅值谱有点不太一样，但是它也可以表达出幅值谱的基频位置（即出现峰值处频率）

octopussheng

我从她的程序里面看就是一个赋值的功率谱计算啊，呵呵，不知道你的想法怎样，也可以帮忙解决一下她的问题哦！

咕噜噜

^_^，我又说错话了？
我没说申老师说得一定正确啊，只是让你们参考一下啊:@L :@L

sssssxxxxx921

功率谱的话  你的程序已经正确了  我再说没什么意义了  呵呵
幅值谱的纵坐标只要去掉功率谱中的平方就可以了   
但是如果纵坐标要表示真正的幅值的话  我也不太清楚这样是不是
但是这样的幅值谱也是正确的，而且功率谱也有很多画法
在绘制功率普的时候最好用matlab自带的命令就可以了   这我也是刚从一个高手那学的  哈哈
幅值谱不知道有没有   
另外  小咕  我是不是有责备你的意思啊 哈哈   大家在一起只是交流， 你对别人的问题非常热心
而且学识也很渊博，很佩服你的     你那样说   反到让我不好意思了:lol

秋月

谢谢啊!
这个程序看不太懂,还请指点power=abs(Y(1:N/2)).^2;
freq=(1:N/2)/(N/2)*nyquist;   ％求频率
为什么要平方?

秋月

采样定理：
当采样频率fs.max大于信号中，最高频率fmax的2倍时，即：fs.max>=2fmax,则采样之后的数字信号完整地保留了原始信号中的信息，一般取2.56-4倍的信号最大频率；采样定理又称奈奎斯特定理。
修改之后是这样的,麻烦帮忙看看理解的对不对?
打扰了,谢谢!
clc
clear all
x0=[0;1.0;1 ];
[t,y]=ode45('duffing',[0:0.01:100],x0);
%频谱图
fs=300;          %采样频率
N=50;           %采样点数
x=y(5000:100:10000,1);  %生成信号（Y轴方向的位移）
Y=fft(x,N);
Pyy = Y.* conj(Y) /N;
df=fs/(N-1);          %频率分辨率Hz
f=(0:N-1)*df;    %频域序列
figure
% plot(f,Pyy);
% figure
plot(f(1:fix(N/2)),Pyy(1:fix(N/2)));
xlabel('Frequency(HZ)');
ylabel('Magnitude');
title('频率幅值图');
hold on

octopussheng

这个意思就是求功率谱啊，可以看功率谱的定义的，呵呵！

秋月

麻烦你简单介绍一下功率谱和幅植谱的区别?
怎么做才是幅植谱?
谢谢!