搜索
 找回密码
 立即注册
查看: 499|回复: 4

【源程序】时域波转换成频域波(傅里叶变换)

[复制链接]
  • TA的每日心情
    无聊
    2017-10-17 10:22
  • 签到天数: 172 天

    [LV.7]常住居民III

    37

    主题

    337

    帖子

    4130

    积分

    管理员

    Rank: 9Rank: 9Rank: 9

    积分
    4130
    发表于 2017-1-13 16:54:21 | 显示全部楼层 |阅读模式
    1. %FFT变换,获得采样数据基本信息,时域图,频域图
    2. %这里的向量都用行向量,假设被测变量是速度,单位为m/s
    3. clear;
    4. close all;

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

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

    23. %显示原始数据曲线图(时域)
    24. subplot(2,1,1);
    25. semilogx(x,y,'k','LineWidth',2)                                                                                %显示原始数据曲线图
    26. axis([min(x) max(x) 1.1*floor(min(y)) 1.1*ceil(max(y))])               %优化坐标,可有可无
    27. xlabel('Time ');
    28. ylabel('Pressure [Pa]');
    29. title('Time domain');
    30. grid on;

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

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

    40. %显示频谱图(频域)
    41. subplot(2,1,2)
    42. stem(f(1:round(N/2)),Mag(1:round(N/2)),'r','LineWidth',2,'Marker','none')                            %显示频谱图,后面是表示去除标杆上的圈
    43. %stem(f(1:round(N/2)),Mag(1:round(N/2)),'r','LineWidth',2)
    44. %                 |
    45. %                将这里的Pyy改成Mag就是 幅值-频率图了
    46. axis([min(f(1:round(N/2))) max(f(1:round(N/2))) 1.1*floor(min(Mag(1:round(N/2)))) 1.1*ceil(max(Mag(1:round(N/2))))])
    47. xlabel('Frequency [Hz]')
    48. ylabel('Pressure Magnitude')
    49. title('Frequency domain')
    50. grid on;
    51. h=figure(1);
    52. axis;
    53. h_axis=get(h,'Children');
    54. set(h_axis,'LineWidth',2); %坐标轴加粗

    55. %返回最大能量对应的频率和周期值
    56. [a b]=max(Mag(1:round(N/2)));
    57. fprintf('\n傅立叶变换结果:\n')
    58. fprintf('           FFT_f = %1.3f Hz\n',f(b))               %输出最大值对应的频率
    59. fprintf('           FFT_T = %1.3f s\n',1/f(b))             %输出最大值对应的周期

    复制代码


  • TA的每日心情
    无聊
    2017-10-17 10:22
  • 签到天数: 172 天

    [LV.7]常住居民III

    37

    主题

    337

    帖子

    4130

    积分

    管理员

    Rank: 9Rank: 9Rank: 9

    积分
    4130
     楼主 发表于 2017-1-13 20:44:59 | 显示全部楼层
    菠菜有料 发表于 2017-1-13 18:40
    用插入代码这个功能,有个好处是可以直接复制代码,

    :victory::victory::victory:
  • TA的每日心情
    难过
    2017-12-22 18:40
  • 签到天数: 111 天

    [LV.6]常住居民II

    17

    主题

    173

    帖子

    1330

    积分

    管理员

    学霸

    Rank: 9Rank: 9Rank: 9

    积分
    1330
    发表于 2017-1-13 18:40:46 | 显示全部楼层
    用插入代码这个功能,有个好处是可以直接复制代码,
    [fly]我爱学习,哈哈[/fly]
  • TA的每日心情
    难过
    昨天 15:59
  • 签到天数: 395 天

    [LV.9]以坛为家II

    113

    主题

    1071

    帖子

    1万

    积分

    版主

    Rank: 7Rank: 7Rank: 7

    积分
    14314

    荣誉会员学者

    QQ
    发表于 2017-1-13 22:35:17 | 显示全部楼层
    感谢分享!
    流场,温度场,低温等离子体
    回复

    使用道具 举报

  • TA的每日心情

    2017-12-3 10:58
  • 签到天数: 1 天

    [LV.1]初来乍到

    0

    主题

    5

    帖子

    29

    积分

    新手上路

    Rank: 1

    积分
    29
    发表于 2017-12-3 11:32:09 | 显示全部楼层
    楼楼,可以的

    发表回复

    您需要登录后才可以回帖 登录 | 立即注册

    本版积分规则

    4