大数据采集与信号处理

数据信息采集与处理

基本内容：基于FFT的功率谱分析程序设计与应用

1.基本要求

1)对一个人为产生的信号进行采用FFT变换方法进行功率谱分析。

已知信号x(n)=80.0*COS(2*3.14*SF*n/FS)

式中： n=0,1,2 ……N-1

SF---信号频率

FS---采样频率

其FFT变换结果X(k)可用下面提供的FFT子程序求出，计算功率谱的公式为： W（k）=2(XR(k)2 +XI(k)2)/N

式中：k=0,1,2 ……N/2-1

XR(k)--- X(k)的实部

XI(k)--- X(k)的虚部

请用VB，VC或C++Builder编译器编程，或采用MATLAB计算，或采用高级语言调用MATLAB计算。处理结果为采用窗口显示时域波形和频域波形。

此信号的时域谱、频域谱、功率谱如下面图1~图3所示：

图1

图2

图3 其MATLAB代码为：

FS=200;

SF=10;

N=1024;

n=0:N-1;

t=n/FS;

x=80.0*cos(2*3.14*SF*t);

figure;

plot(t，x);

xlabel('t');

ylabel('y');

title('x=80.0*cos(2*3.14*SF*t)时域波形');

grid;

y=fft(x,N);

mag=abs(y);

f=(0:length(y)-1)*FS/length(y);%进行对应的频率转换

figure;

plot(f(1:N/2),mag(1:N/2));%做频谱图

xlabel('频率(Hz)');

ylabel('幅值');

title('x=80.0*cos(2*3.14*SF*t)幅频谱图N=1024');

grid;

Py =2*(y.*conj(y))/N; %计算功率谱密度Py

figure;

plot(f(1:N/2),Py(1:N/2));

xlabel('频率(Hz)');

ylabel('功率谱密度');

title('x=80.0*cos(2*3.14*sf*t)功率谱密度');

grid;

2)对实验所采集的转子振动信号进行频谱分析

图4

其MATLAB代码：

SF=1000;

fid = fopen('D:\数据采集与信号处理作业\信号处理大作业

\sanjiao_45HZ_1024_1000HZ');%转子信号

[a,N]= fscanf(fid,'%f');

fclose(fid);

y=fft(a,N);%FFT运算

Pyy =sqrt(y.*conj(y))*2.0/N; %取功率普密度

f=(0:length(Pyy)-1)*SF/length(Pyy);

LPyy=20*log10(Pyy);

plot(f(1:N/2),Pyy(1:N/2));%输出FS/2点幅频谱图

xlabel('频率(Hz)');

ylabel('幅值');

title('转子振动信号频谱图');

grid;

2.讨论

1)信号经过均值化处理或不经过均值化处理的结果比较：

图5

图6

通过以上两个图（图5、图6）的分析，我们可以看出均值化处理后的频谱的低频段消失，这就去去除了常规的干扰频谱，如环境噪声等，对我们进行频谱分析有很大作用。

其MATLAB代码为：

Fs=200; %采样频率

n=0:1:200;

N=1024;

SF=10; %信号频率

xn=80.0*cos(2*3.14*SF*n/Fs)+70;%产生波形序列

window=boxcar(length(xn)); %矩形窗

nfft=512;%采样点数

[Pxx,f]=periodogram(xn,window,nfft,Fs); %直接法

plot(f,Pxx);

xlabel('频率(Hz)');

ylabel('功率谱密度');

title('含有直流分量的余弦曲线未均值化的功率谱波形图');

grid;

2)采用不同窗函数时的谱结果（矩形窗函数, 汉宁窗函数，汉明窗）如图7~12所示。

图7

图8

图9

图10

图11

图12

其MATLAB代码为：

y=fft(x,N);

mag=abs(y);

f=(0:length(y)-1)*FS/length(y);%进行对应的频率转换

w_han=(hanning(N))';

y1=x.*w_han;

figure;

plot(t,y1);

xlabel('t');

ylabel('y');

title('汉宁窗时域波形');

grid;

y2=mag.*w_han;

figure;

plot(f(1:N/2),y2(1:N/2));

xlabel('频率(Hz)');

ylabel('幅值');

title('汉宁窗频域特性');

grid;

w_rect=(rectwin(N))';

y3=x.*w_rect;

figure;

plot(t,y3);

xlabel('t');

ylabel('y');

title('矩形窗时域波形');

grid;

y4=mag.*w_rect;

figure;

plot(f(1:N/2),y4(1:N/2));

xlabel('频率(Hz)');

ylabel('幅值');

title('矩形窗频域特性');

grid;

w_ham=(hamming(N))';

y5=x.*w_ham;

figure;

plot(t,y5);

xlabel('t');

ylabel('y');

title('汉明窗时域波形');

grid;

y6=mag.*w_ham;

figure;

plot(f(1:N/2),y6(1:N/2));

xlabel('频率(Hz)');

ylabel('幅值');

title('汉明窗频域特性');

grid;

3)典型函数的频谱（矩形窗函数, 汉宁窗函数，直线，阶跃函数，δ函数，方波，三角波等），如图13~18所示。

图13

图 14

图

15

图 17

图18 此部分MATLAB代码如下：

t=0:0.001:0.2;

N=256;

FS=300;

w=boxcar(N); %产生信号

figure;

plot(w);title('矩形窗函数的时域波形图');

axis([0,260,0,2]);grid on;

y=fft(w,N); %FFT运算

mag=abs(y);%取幅值

f=(0:length(y)-1)*FS/length(y);

figure;

plot(f(1:N/2),mag(1:N/2)); %输出FS/2点幅频谱图title('矩形窗函数频域波形图');grid;

xlabel('频率');ylabel('幅值');

t=0:0.001:0.2;

N=256;

FS=300;

w=hanning(N); %产生信号

figure;

plot(w);title('汉宁窗函数的时域波形图');

grid on;

y=fft(w,N); %FFT运算

mag=abs(y); %取幅值

f=(0:length(y)-1)*FS/length(y);

figure;

plot(f(1:N/2),mag(1:N/2)); %输出FS/2点幅频谱图title('汉宁窗函数频域波形图');

grid on;

xlabel('频率');ylabel('幅值');

t=0:0.001:0.2;

N=256;

FS=300;

w=1; %产生信号

y=fft(w,N); %FFT运算

mag=abs(y); %取幅值

f=(0:length(y)-1)*FS/length(y);

figure;

plot(f(1:N/2),mag(1:N/2)); %输出FS/2点幅频谱图

title('直线频域波形图');

grid on;

xlabel('Frequency(Hz)');ylabel('Magnitude');%阶跃函数的频域波图

clc;clf;t=0:0.001:0.2;

N=256;

FS=300;

w=ones(1,N); %产生信号

figure;

plot(w);

title('阶跃函数的时域波形图');

grid on;

y=fft(w,N); %FFT运算

mag=abs(y);%取幅值

f=(0:length(y)-1)*FS/length(y);

figure;

plot(f(1:N/2),mag(1:N/2)); %输出FS/2点幅频谱图

title('阶跃函数的频域波形图');

grid on;

xlabel('频率');ylabel('幅值');

t=0:0.001:0.2;

N=256;

FS=300;

w=zeros(1,N);w(1)=1; %产生信号

figure;

plot(w);

grid on;

title('δ函数的时域波形图');

y=fft(w,N);%FFT运算

mag=abs(y);%取幅值

f=(0:length(y)-1)*FS/length(y);

figure;

plot(f(1:N/2),mag(1:N/2)); %输出FS/2点幅频谱图

title('δ函数的频域波形图');

grid on;

xlabel('频率');ylabel('幅值');

t=0:0.001:0.2;

N=256;

FS=300;

w=square(2*pi*50*t); %产生信号

figure;

plot(t,w);

title('方波的时域波形图');

axis([0,0.2,-0.2,1.2]);

grid on;

y=fft(w,N); %FFT运算

mag=abs(y); %取幅值

f=(0:length(y)-1)*FS/length(y);

figure;

plot(f(1:N/2),mag(1:N/2)); %输出FS/2点幅频谱图title('方波的频域波形图');

grid on;

xlabel('频率');ylabel('幅值');

t=0:0.001:0.2;

N=256;

FS=300;

w=sawtooth(2*pi*50*t,0.5);

figure;

plot(t,w);

grid on;

title('三角波的时域波形图');%产生信号

y=fft(w,N); %FFT运算

mag=abs(y); %取幅值

f=(0:length(y)-1)*FS/length(y);

figure;

plot(f(1:N/2),mag(1:N/2)); %输出FS/2点幅频谱图title('三角波的频域波形图');

grid on;

xlabel('频率');ylabel('幅值');

4)整周期和非整周期采样时两者的比较，如图19~图22所示。

图 19

图 20

图 21

图 22 此部分的MATLAB代码为：

FS=900;

SF=10;

N=1024;

n=0:N-1;

t=n/FS;

x=80.0*cos(2*3.14*SF*t);

figure;

plot(t,x);

xlabel('t');

ylabel('y');

title('非整周期时域波形');

grid;

y=fft(x,N);

mag=abs(y);

f=(0:length(y)-1)*FS/length(y);%进行对应的频率转换figure;

plot(f(1:N/2),mag(1:N/2));%做频谱图

xlabel('频率(Hz)');

ylabel('幅值');

title('非整周期幅频谱图N=1024');

grid;

FS=1024;

SF=10;

N=1024;

n=0:N-1;

t=n/FS;

x=80.0*cos(2*3.14*SF*t);

figure;

plot(t,x);

xlabel('t');

ylabel('y');

title('整周期时域波形');

grid;

y=fft(x,N);

mag=abs(y);

f=(0:length(y)-1)*FS/length(y);%进行对应的频率转换figure;

plot(f(1:N/2),mag(1:N/2));%做频谱图

xlabel('频率(Hz)');

ylabel('幅值');