数据采集与信号处理(个人整理版)

一．基本内容：基于FFT的功率谱分析程序设计与应用
1）对一个人为产生的信号，采用FFT变换方法进行功率谱分析
已知信号x(n)=100.0*COS(2*3.14*S F*n/F S) +50.0*COS(4*3.14*S F*n/F S)，n=0,1,2 ……N-1
式中：N---采样点数，必须为2x
S F---信号频率，S F=10 Hz
F S---采样频率
其FFT变换结果X(k)可用下面提供的FFT子程序求出，计算功率谱的公式为：
G(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计算。

处理结果（时域波形或频谱）采用窗口显示，并且拷贝至作业稿中。

注意：除FFT分析程序可参考附件所附的程序外，其余程序必须自己设计（小波或EMD除外），后面限选题目要求与此相同。

该信号的处理结果如下：
其matlab程序如下：
SF=10.0;%信号频率
N=1024;%采样点数
n=0:1023;%时间轴N个点
FS=300;%采样频率
t=n/FS;
x=100.0*cos(2*3.14*SF*t) +50.0*cos(4*3.14*SF*t);
figure;
plot(t,x);%信号输出
xlabel('t');
ylabel('y');
title('时域波形');
grid;
y1=fft(x,N);%FFT运算
mag=abs(y1);
f=(0:length(y1)-1)*FS/length(y1);
figure;
plot(f,mag);%做频谱图
xlabel('频率(Hz)');
ylabel('幅值');
title('幅频谱图');
grid;
Pyy1 =(y1.*conj(y1))*2/N;%取功率普密度
figure;
plot(f,Pyy1);
xlabel('频率(Hz)');
ylabel('功率谱密度');
title('功率谱密度');
grid;
2）对一个用A/D数据采集板采集的信号进行频谱分析参考数据文件名(见作业文件压缩包)：
①fanbo_45HZ_1024_1000HZ，45Hz方波发生器输出信号采集结果，采样频率1000Hz，采样点数1024
② sin_45HZ_1024_1000HZ，45Hz正弦波发生器输出信号采集结果，采样频率1000Hz，采样点数1024
③ A_1_1024_1000HZ，汽轮机轴向振动位移信号采集结果，采样频率1000Hz，采样点数1024
数据文件名的含义，如文件名为fanbo_45HZ_1024_1000HZ，表示45Hz方波信号，采样点数1024，采样频率1000Hz
①该信号处理结果如下：
其matlab程序如下：
SF=45;%信号频率
N=1024;%采样点数
n=0:1023;%时间轴N个点
FS=1000;%采样频率
t=n/FS;
x=load('D:\fanbo_45HZ_1024_1000HZ.txt');
figure;
plot(t,x);%信号输出
xlabel('t');
ylabel('y');
title('时域波形');
grid;
y1=fft(x,N);%FFT运算
mag=abs(y1);
f=(0:length(y1)-1)*FS/length(y1); figure;
plot(f,mag);%做频谱图
xlabel('频率(Hz)');
ylabel('幅值');
title('幅频谱图');
grid;
Pyy1 =(y1.*conj(y1))*2/N;%取功率普密度figure;
plot(f,Pyy1);
xlabel('频率(Hz)');
ylabel('功率谱密度');
title('功率谱密度');
grid;
②该信号的处理结果如下：
其matlab程序基本同上。

③该信号处理结果如下：
其matlab程序基本同上。

3）讨论
1. 信号经过零均值化处理或不经过零均值化处理的结果比较。

对于已知信号x=100.0*cos(2*3.14*SF*t) +50.0进行处理，结果如下：
该信号的非零均值相当于在此信号的低频段叠加上了一个直流分量，可以看到在零频率处出现一个很大的谱峰，并会影响在零频率左、右处的频谱曲线，使之产生较大误差。

零均值化处理后零频率出的谱峰消失，并去除了干扰。

其matlab处理程序如下：
FS=300;%采样频率
n=0:300;
N=1024;%采样点数
SF=10.0;%信号频率
t=n/FS;
x=100.0*cos(2*3.14*SF*t) +50.0;%产生波形序列
window=boxcar(length(x));%矩形窗
[P,f]=periodogram(x,window,N,FS);%直接法
subplot(211)
plot(f,P);
title('信号未零均值化的功率谱波形图');
xlabel('频率');
ylabel('功率谱密度');
%零均值化的信号
avg=mean(x);
x1=x-avg;
window=boxcar(length(x));%矩形窗
[P,f]=periodogram(x1,window,N,FS);%直接法
subplot(212);
plot(f,P);
title('信号零均值化后的功率谱波形图');
xlabel('频率');
ylabel('功率谱密度');
2．采用不同窗函数时的谱结果分析（矩形窗函数, 汉宁窗函数，汉明窗）。

加矩形窗的功率谱图：
加汉宁窗的功率谱图：
加汉明窗的功率谱图：
其matlab程序如下：
FS=300;%采样频率
SF=10;
N=1024;%采样点数
t=0:1/FS:1;
x=100.0*cos(2*3.14*SF*t) +50.0*cos(4*3.14*SF*t);%产生
信号序列
window1=boxcar(length(x));%矩形窗函数
window2=hanning(length(x));%汉宁窗
window3=hamming(length(x));%汉明窗
noverlap=0;
p=0.9;
[P1,f1]=SPECTRUM(x,N,FS,window1,noverlap,p);
index=0:round(N/2-1);
k=index*FS/N;
plot_P1=P1(index+1);
figure;
plot(k,plot_P1);
grid;
title('矩形窗');
[P2,f2]=SPECTRUM(x,N,FS,window2,noverlap,p);%计算序
列的PSD
plot_P2=P2(index+1);
figure;
plot(k,plot_P2);
grid;
title('汉宁窗');
[P3,f3]=SPECTRUM(x,N,FS,window3,noverlap,p);%计算序
列的PSD
plot_P3=P3(index+1);
figure;
plot(k,plot_P3);
grid;
title('汉明窗');
3．人为产生典型信号并进行谱分析，自由讨论计算结果（矩形窗函数, 汉宁窗函数，直线，δ函数，方波等）。

矩形窗函数的时域频域波形图：
汉宁窗函数的时域频域波形图：
直线的时域频域波形图：
δ函数的时域频域波形图：
方波的时域频域波形图：
其matlab程序如下：
%矩形窗函数的频域波形图
N=256;
FS=300;
w=boxcar(N);%产生信号
figure(1);
subplot(211);
plot(w);
title('矩形窗函数的时域波形图');
axis([0,260,0,2]);
grid;
y=fft(w,N);%FFT运算
mag=abs(y);%取幅值
f=(0:length(y)-1)*FS/length(y); subplot(212);
plot(f(1:N/2),mag(1:N/2));
title('矩形窗函数频域波形图'); grid;
xlabel('频率');
ylabel('幅值');
%汉宁窗函数的频域波形图
N=256;
FS=300;
w=hanning(N);%产生信号
figure(2);
subplot(211);
plot(w);
title('汉宁窗函数时域波形图'); grid;
y=fft(w,N);%FFT运算
mag=abs(y);%取幅值
f=(0:length(y)-1)*FS/length(y); subplot(212);
plot(f(1:N/2),mag(1:N/2));
title('汉宁窗函数频域波形图'); grid;
xlabel('频率');
ylabel('幅值');
%直线的频域波形图
N=256;
FS=300;
w=1;%产生信号
figure(3);
subplot(211);
plot(w);
title('直线函数时域波形图'); grid;
y=fft(w,N);%FFT运算
mag=abs(y);%取幅值
f=(0:length(y)-1)*FS/length(y); subplot(212);
plot(f(1:N/2),mag(1:N/2));
title('直线函数频域波形图'); grid;
xlabel('频率');
ylabel('幅值');
%δ函数的频域波形图
N=256;
FS=300;
w=zeros(1,N);w(1)=1;%产生信号figure(4);
subplot(211);
plot(w);
title('δ函数时域波形图'); grid;
y=fft(w,N);%FFT运算
mag=abs(y);%取幅值
f=(0:length(y)-1)*FS/length(y); subplot(212);
plot(f(1:N/2),mag(1:N/2));
title('δ函数频域波形图'); grid;
xlabel('频率');
ylabel('幅值');
%方波的频域波形图
t=0:0.001:0.2;
N=256;
FS=300;
w=square(2*3.14*50*t);%产生信号figure(5);
subplot(211);
plot(t,w);
title('方波的时域波形图');
axis([0,0.2,-0.2,1.2]);
grid;
y=fft(w,N);%FFT运算
mag=abs(y);%取幅值
f=(0:length(y)-1)*FS/length(y); subplot(212);
plot(f(1:N/2),mag(1:N/2));
title('方波的频域波形图'); grid;
xlabel('频率');
ylabel('幅值');
4．整周期和非整周期采样时两者频谱的比较。

其matlab程序如下：
%汉明窗周期采样
FS=300;%采样频率
n=0:1023;
N=1024;%采样点数
SF=20;%信号频率
t=n/FS;
x=100.0*cos(2*3.14*SF*t) +50.0*cos(4*3.14*SF*t);
window=hamming(length(x));%汉明窗
[P1,f1]=periodogram(x,window,N,FS);%直接法
subplot(211);
plot(f1,P1);
grid;
title('信号周期采样');
xlabel('频率');
ylabel('幅值');
%汉明窗非周期采样
T=1.6*n;
x1=0.5*(1-cos(2*3.14*T/N));
z=x.*x1;
y=fft(z,N);
mag=abs(y);
f=(0:length(y)-1*FS/length(y));
subplot(212);
plot(f(1:N/2),mag(1:N/2));
xlabel('频率');
ylabel('幅值');
title('信号非周期采样');
grid;
二、自选内容
1) 设计计算包络谱的程序，并对一个模拟幅值调制信号进行包络谱分析
信号可自行编程产生，信号为：
12()120.0cos(2)(10.85sin(2))x t f t f t ππ=+
其中：f 1=100 Hz ，f 2=10 Hz 。

要求：包络谱计算过程不可调用现成的函数（如MATLAB 相关函数）。

其matlab程序如下：
FS=300;
SF1=100.0;
SF2=10.0;
N=1024;
n=0:1023;
t=n/FS;
x=120.0*cos(2*3.14*SF1*t)*(1+0.85*sin(2*3.14*SF2*t))
';
y=fft(x,N);
mag=abs(y);
f=(0:length(y)-1)*FS/length(y);%进行对应的频率转换w1=(hanning(N))';
y1=x.*w1;
subplot(311);
plot(t,y1);
xlabel('t');
ylabel('y');
title('汉宁窗时域波形');
grid;
y11=mag.*w1;
subplot(312);
plot(f(1:N/2),y11(1:N/2));
xlabel('频率');
ylabel('幅值');
title('汉宁窗频域特性');
grid;
y2=hilbert(y1);% hilbert变换
y3=abs(y2);% 绝对值
y=y3-mean(y3);%均值化处理
p=abs(fft(y3,N));
f=(0:length(fft(y3,N))-1)*FS/length(fft(y3,N)); subplot(313);
plot((0:N/2-1)/N*FS,p(1:N/2));
xlabel('频率');
ylabel('包络谱幅值');
title('汉宁窗包络谱频域特性');
grid;。