频域分析一维离散信号

第二种方法频域分析

clear all; close all;
clc;
%FFT变换,获得采样数据基本信息，时域图，频域图
%这里的向量都用行向量，假设被测变量是速度，单位为m/s
clear;
close all;
data=importdata(‘TEK0001.csv’); %读取csv数据文件
% disp(data0); %disp函数：显示文本或数组
循环读取离散信号数值，后来觉得不适用于千数量级的数据读取
% for i=1:2250
% v(i,1)=data(i,4);
% figure(1)
% plot(voltagechange)
% xlabel(‘数目’),ylabel(‘y/voltagechange’); %添加标签
% hold on;
% end
数据的行列变换
A=data; %将测量数据赋给A，此时A为N×2的数组
x=A(:,1); %将A中的第一列赋值给x，形成时间序列
x=x’; %将列向量变成行向量
y=A(:,4); %将A中的第二列赋值给y，形成被测量序列
y=y’; %将列向量变成行向量

% m=15;%表示平滑滤波窗长度,这是长度为奇数的情况
% %前m/2,最后m/2个点没滤波，设为原来的值就行
% for i=1:length(y)-m+1
% y(i+(m-1)/2)=sum(y(i:i+m-1))/m;
% end
% plot(y);
% title(‘原始的信号’)
% figure(2)
% [thr,sorh,keepapp]=ddencmp(‘den’,‘wv’,y);
% xd=wdencmp(‘gbl’,y,‘db4’,2,thr,sorh,keepapp);
% xlabel(‘时间 t/s’);ylabel(‘频率 f/Hz’);
% plot(xd);
% title(‘消除噪声的信号’)

%显示数据基本信息
fprintf(’\n数据基本信息：\n’)
fprintf(’ 采样点数 = %7.0f \n’,length(x)) %输出采样数据个数
fprintf(’ 采样时间 = %7.3f s\n’,max(x)-min(x)) %输出采样耗时
fprintf(’ 采样频率 = %7.1f Hz\n’,length(x)/(max(x)-min(x))) %输出采样频率
fprintf(’ 最小速度 = %7.3f m/s\n’,min(y)) %输出本次采样被测量最小值
fprintf(’ 平均速度 = %7.3f m/s\n’,mean(y)) %输出本次采样被测量平均值
fprintf(’ 速度中值 = %7.3f m/s\n’,median(y)) %输出本次采样被测量中值
fprintf(’ 最大速度 = %7.3f m/s\n’,max(y)) %输出本次采样被测量最大值
fprintf(’ 标准方差 = %7.3f \n’,std(y)) %输出本次采样数据标准差
fprintf(’ 协 方 差 = %7.3f \n’,cov(y)) %输出本次采样数据协方差
fprintf(’ 自相关系数 = %7.3f \n\n’,corrcoef(y)) %输出本次采样数据自相关系数

%傅立叶变换
y=y-mean(y); %消去直流分量，使频谱更能体现有效信息
Fs= 2501.1; %得到原始数据data.txt时，仪器的采样频率。就是length(x)/(max(x)-min(x));
N=2250; %data.txt中的被测量个数，即采样个数。其实就是length(y);
z=fft(y);

%频谱分析
f=(0:N-1)Fs/N;
Mag=2abs(z)/N; %幅值，单位同被测变量y
Pyy=Mag.^2; %能量；对实数系列X，有 X.*X=X.*conj(X)=abs(X).^2=X.2，故这里有很多表达方式

%显示原始数据曲线图（时域）
subplot(2,1,1);
plot(y) %显示原始数据曲线图
%axis([min(x) max(x) 1.1floor(min(y)) 1.1ceil(max(y))]) %优化坐标，可有可无
xlabel(‘时间 (s)’);
ylabel(‘被测变量y’);
title(‘处理后信号(时域)’);
grid on;
%显示频谱图(频域)
subplot(2,1,2)
plot(f(1:N/2),Mag(1:N/2),‘r’) %显示频谱图 |
% 将这里的Pyy改成Mag就是 幅值-频率图了
% axis([min(f(1:N/2)) max(f(1:N/2)) 1.1*floor(min(Pyy(1:N/2))) 1.1 ceil(max(Pyy(1:N/2)))])
xlabel(‘频率 (Hz)’)
ylabel(‘能量’)
title(‘频谱图(频域)’)
grid on;

%返回最大能量对应的频率和周期值
[a, b]=max(Pyy(1:N/2));
fprintf(’\n傅立叶变换结果：\n’)
fprintf(‘FFT_f = %1.3f Hz\n’,f(b)) %输出最大值对应的频率
fprintf(‘FFT_T = %1.3f s\n’,1/f(b)) %输出最大值对应的周期

Hd =lowpass_ly;%此处是自己设计的一个低通滤波器通过（filterDesigner）
LPF_Data = filter(Hd,y);
figure;plot(LPF_Data)
title(‘低通滤波之后的波形’)
%-------------------低通滤波之后频谱分析
SampleFre=2250;
FFT_LPF_Data = fft(LPF_Data);
Amplitude_LPF = abs(FFT_LPF_Data);
Amplitude_LPF = 2Amplitude_LPF/length(Amplitude_LPF);
% Amplitude_LPF(2:end) = 2Amplitude_LPF(2:end);
Frequence = (0:(length(Amplitude_LPF)/2-1))/length(Amplitude_LPF)*SampleFre;
figure;plot(Frequence,Amplitude_LPF(1:length(Frequence)))
title(‘低通滤波之后的频谱’)
xlabel(‘t(s)’)
ylabel(‘Voltage(v)’)

滤波结果存在两个问题，fir线性滤波存在前端数据异常情况，再就是参数需要调整的，调试起来也挺麻烦的，不过图片上滤的有点过了，这回学会了用filterDesigner设计滤波器，其实原理并不难。
是filter的问题，filter是将两个序列进行卷积之后取前面那一段。比如序列x1长度为L1,序列x2长度为L2,那么将两个序列卷积之后的长度为L1+L2-1,如果x1为信号，x2为滤波器，那么系统输出是需要L1即信号长度的输出序列，filter自动截取的是前L1个序列，将filter函数改成两个函数卷积conv之后取中间L1长度输出，应该就不会有延迟了。