开篇语
大学嘛~主业就念书考试做实验对头不对头?念书这个不用说,考试这玩意也是家常便饭了,唯一有点搞头的也就做实验了,然后我们学院大三的,经历过考完机械原理的快乐不到一天以后,就迎来了工程测试技术带来的绝望!!中午去打印店,黑漆漆一片都是在打印频谱图的,盛况空前٩꒰´`꒱۶ෆ͙⃛
正文
食用本文前请先细读上图文字,有利于前后文发展~~~
出现镇文神图的现象~究其原因,是因为,在机械原理考试之后,有一门神奇的课程实验,不,更准确的说是有一个神奇的老师,给我们留下了一个很蛋疼的实验报告。所以,我们疯狂了,寝室里面最爱打游戏的人都活生生的学到两点才意犹未尽的停下了写Matlab,准备第二天早早的起来写程序,所以,历经千辛万苦,我觉得我还要把辛辛苦苦写出来的源码放到简书上,哪怕为了纪念一番那也是好的!!
源码
谁说只有软件项目能开源?不知道以前有没有人把作业开源的,但是对于这个磨人的实验,我今天要把它放到简书上,说不定我们的下一届,能在未来上简书搜索互相关,Matlab 的时候能找到他们的前辈留下来的这份报告,这是我们曾经与黄教主(实验老师)战斗留下的可歌可泣的传说!!
废话不说的,代码先行:
fs=400; %定义采样频率
N=256; %定义采样点数
n=0:N-1;t=n/fs; %在采样时间内按时间均分地描绘128个点
x1=sin(10*pi*t)+2*sin(20*pi*t)+3*sin(30*pi*t); %信号产生函数
plot(t,x1);title('原始信号'),xlabel('时间/s');hold; %描绘出原始信号曲线
y=fft(x1,N); %傅里叶离散变换
y1=y/256;
y2=y1*2;
mag=abs(y2); %求y2的模
f=(0:128)*400/256;
plot(f,mag(1:129));title('频域信号'),xlabel('时间/s');
py1=angle(y);
py2=180/pi*py1;
plot(f,py1(1:129));title('相位'),xlabel('频率/f') %绘出相频谱
plot(f,mag(1:129));title('幅值'),xlabel('频率/f') %绘出幅频谱
sp=real(y);
plot(f,sp(1:129));title('实频谱'),xlabel('频率/f') %绘出实频谱
xp=imag(y);
plot(f,xp(1:129));title('虚频谱'),xlabel('频率/f') %绘出虚频谱
power=mag(1:129).^2;
plot(f,power);title('功率谱'),xlabel('频率/f') %绘制功率谱
lag=200;
[c,lags]=xcorr(x1,lag,'unbiased');
plot(lags/fs,c);title('自相关信号'),xlabel('时间/s')
lag=200;
[c,lags]=xcorr(x1,lag,'unbiased');
m=length(c)-1;
z=fft(c,m);
mag_z=abs(z);
ff=(0:m-1)*fs/m;
plot(ff(1:m/2),mag_z(1:m/2)*2/m);title('自相关频谱'),xlabel('频率/Hz');
lag=200;
x2=wgn(1,256,1)
[lagc1,lags1]=xcorr(x1,x2,lag,'unbiased');
m1=length(lagc1)-1;
z1=fft(lagc1,m1);
mag_z1=abs(z1);
ff1=(0:m1-1)*fs/m1;
plot(lags1/fs,lagc1);title('互相关信号'),xlabel('时间/s');
plot(ff1(1:m1/2),mag_z1(1:m1/2)*2/m1);title('互相关频谱'),xlabel('频率/Hz');
后来的学妹们,这一份是草稿,也许其中会有参数定义错误,因为这一份代码有四个学长用过,每个人都留下了自己的痕迹,当然,主要还是我找的那一份不知何年何月在此作战的学长留下的宝贵代码~~
下面是我的手稿:这一份亲测可行,我还会放上我的实验图。
实验数据图片
结束语
致我可爱的学弟学妹们,如果你们实在走投无路,不知如何下手,可用我的实验报告以之为参考,但是学习是你自己的,我还是希望你们能够独立完成,如果时间充裕,情况允许,你还是得自己写,Matlab是一个强大的工具,学会这个东西对你好处非凡,加油!