zoukankan      html  css  js  c++  java
  • 离散信号MATLAB频谱分析程序

    from http://blog.csdn.net/u012129372/article/details/26565611

    %FFT变换,获得采样数据基本信息,时域图,频域图
    %这里的向量都用行向量,假设被测变量是速度,单位为m/s
    clear;
    close all;

    load data.txt              %通过仪器测量的原始数据,存储为data.txt中,附件中有一个模版(该信号极不规则)
    A=data;                                        %将测量数据赋给A,此时A为N×2的数组
    x=A(:,1);                                     %将A中的第一列赋值给x,形成时间序列
    x=x';                                           %将列向量变成行向量
    y=A(:,2);                                     %将A中的第二列赋值给y,形成被测量序列
    y=y';                                           %将列向量变成行向量

    %显示数据基本信息
    fprintf(' 数据基本信息: ') 
    fprintf('        采样点数 = %7.0f ',length(x))                         %输出采样数据个数
    fprintf('        采样时间 = %7.3f s ',max(x)-min(x))                    %输出采样耗时
    fprintf('        采样频率 = %7.1f Hz ',length(x)/(max(x)-min(x)))   %输出采样频率
    fprintf('        最小速度 = %7.3f m/s ',min(y))                         %输出本次采样被测量最小值
    fprintf('        平均速度 = %7.3f m/s ',mean(y))                      %输出本次采样被测量平均值
    fprintf('        速度中值 = %7.3f m/s ',median(y))                   %输出本次采样被测量中值
    fprintf('        最大速度 = %7.3f m/s ',max(y))                          %输出本次采样被测量最大值
    fprintf('        标准方差 = %7.3f ',std(y))                               %输出本次采样数据标准差
    fprintf('       协 方 差 = %7.3f ',cov(y))                                %输出本次采样数据协方差
    fprintf('     自相关系数 = %7.3f ',corrcoef(y))                       %输出本次采样数据自相关系数
      
    %显示原始数据曲线图(时域)
    subplot(2,1,1);
    plot(x,y)                                                                                %显示原始数据曲线图
    axis([min(x) max(x) 1.1*floor(min(y)) 1.1*ceil(max(y))])             %优化坐标,可有可无
    xlabel('时间 (s)');
    ylabel('被测变量y');
    title('原始信号(时域)');
    grid on;

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

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

    %显示频谱图(频域)
    subplot(2,1,2)
    plot(f(1:N/2),Pyy(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(' 傅立叶变换结果: ') 
    fprintf('           FFT_f = %1.3f Hz ',f(b))             %输出最大值对应的频率
    fprintf('           FFT_T = %1.3f s ',1/f(b))          %输出最大值对应的周期

    附件data.txt下载地址http://u.115.com/file/e6cqg126

  • 相关阅读:
    nullnullConnecting with WiFi Direct 与WiFi直接连接
    nullnullUsing WiFi Direct for Service Discovery 直接使用WiFi服务发现
    nullnullSetting Up the Loader 设置装载机
    nullnullDefining and Launching the Query 定义和启动查询
    nullnullHandling the Results 处理结果
    装置输出喷泉装置(贪心问题)
    数据状态什么是事务?
    停止方法iOS CGD 任务开始与结束
    盘文件云存储——金山快盘
    函数标识符解决jQuery与其他库冲突的方法
  • 原文地址:https://www.cnblogs.com/gisalameda/p/5530973.html
Copyright © 2011-2022 走看看