matlab中一数组的时域曲线如何转换为频域曲线

Posted

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了matlab中一数组的时域曲线如何转换为频域曲线相关的知识,希望对你有一定的参考价值。

我有一组时间与数值的数组,请问怎么把时间改为频率
我知道用傅里叶变换 但是得出的结果不正确,请大神们会的详细点说

一般情况下得到的离散的数据点没有明显的拟合函数,这时候可以考虑用离散傅里叶变换。matlab中的fft函数可以完成这个功能。

特殊情况下可以看出数据点所满足的解析式,使用拟合,然后对拟合得到的函数进行傅里叶变换,用matlab的fourier函数即可。


由于一般情况中的fft更具有应用性,下面着重举例说明fft。

引用一段matlab帮助文件提供的代码作说明:(%后面是中文或英文注释)

clc;clear;

Fs = 1000;                    % Sampling frequency,取样频率

T = 1/Fs;                      % Sample time,采样时间间隔

L = 1000;                     % Length of signal,总时间

t = (0:L-1)*T;                % Time vector,时间向量

% Sum of a 50 Hz sinusoid and a 120 Hz sinusoid 信号函数,提供50Hz和120Hz的主频率

x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t); 

y = x + 2*randn(size(t));     % Sinusoids plus noise 信号函数加上模拟的噪音

plot(Fs*t(1:100),y(1:100))    % 信号图

title('Signal Corrupted with Zero-Mean Random Noise')

xlabel('time (milliseconds)')

NFFT = 2^nextpow2(L); % Next power of 2 from length of y 扩充采样点,由1000变为1024

Y = fft(y,NFFT)/L;%除一个L,使归一化,可以不除,不影响对主频率的判断

f = Fs/2*linspace(0,1,NFFT/2+1);%这里除以2是因为fft的对称性,因此只画一半

% Plot single-sided amplitude spectrum.

figure

plot(f,2*abs(Y(1:NFFT/2+1))) %乘2是为了归一化,因为右边一半的fft图像没画;不乘,不影响对主频率的判断

title('Single-Sided Amplitude Spectrum of y(t)')

xlabel('Frequency (Hz)')

ylabel('|Y(f)|')


这行代码“f = Fs/2*linspace(0,1,NFFT/2+1);”如果不理解可以写成

“f = Fs*linspace(0,1,NFFT);”然后把后面的plot行的乘2去掉,NFFT/2+1也改成NFFT,这就等于没有折叠的状态。

至于为什么对称、为什么表达式是这样,就需要去做DFT数学推导了,这里不做推导。


得图如下:

折叠了的fft图。在50与120Hz处有明显的主峰。


未折叠的fft图。右边两个峰值并没有实际意义,只是由于对称性而得到的。

参考技术A meshgrid用来生成网格矩阵,简单地讲,就是把给定的x和y中元素的两两组合都生成出来,这样每一对(x,y)再计算一个对应的z,显然这样得到的是一个z的曲面。但该语句不是必须的,有时候我们只想获得一条三维曲线而已,并不想知道所有x, y元素两两组合的结果是什么,组合我们已经定义好了

参考代码:
clc
clear all
close all
tic
n = 10;
x = 1:n; % x坐标
y = 1:n; % y坐标
%%
% meshgrid演示
[X, Y] = meshgrid(x, y); % meshgrid 函数用来生成网格矩阵
Z = X.^2 + Y;
figure
mesh(X, Y, Z);
grid on
xlabel('x');
ylabel('y');
zlabel('z');
%%
% 不用meshgrid的情况
z = x.^2 + y;
figure
% mesh(x, y, z); % 没有meshgrid生成底面矩阵时,该句出错
plot3(x, y, z); % 一组(x, y)对应一个z值,因此x和y元素个数必须一致
grid on
xlabel('x');
ylabel('y');
zlabel('z');
参考技术B 可以使用快速傅里叶变换:fft。
在matlab命令窗口中输入:help fft.
参考技术C fft函数,傅里叶变换 参考技术D 傅立叶变换

第二章:整车发动机激励--快速傅里叶变换

摘要:第一章介绍了缸压曲线和曲轴转速转化成曲轴中心点的时间-集中力;但是这个集中力并不是我们可以直接使用的频域-集中力;需要经过快速傅里叶变换,将时域力转化成频域力。

  • 快速傅里叶变换的理解

    还是老规矩,我们需要什么?我们需要从时域到频域的转换。怎么理解呢?当然是先找一个已知的频率来看看是什么样的。从网上盗了一张动态图,从图中我们已知各个三角函数的周期(频率),可以做出复杂的曲线。【先看成果,我们的目的是从复杂的曲线中分离出各个周期(频率)】

import numpy as np
from scipy.fftpack import  fft,ifft
import matplotlib.pyplot as plt
import seaborn

"采样点选取,0-2间取2000个采样点,时域横坐标单位s"
x=np.linspace(0,2,1000)
print(x)
"构造信号曲线,其中包含频率1Hz,2Hz,3Hz,4Hz"
y=2*np.sin(2*np.pi*1*x)/np.pi-2*np.sin(2*np.pi*2*x)/(2*np.pi)+2*np.sin(2*np.pi*3*x)/(3*np.pi)\\
  -2*np.sin(2*np.pi*4*x)/(4*np.pi)

plt.plot(x, y)
plt.show()
已知频率信号曲线绘制

 

 

  • 快速傅里叶变换

 

import numpy as np
from scipy.fftpack import  fft,ifft
import matplotlib.pyplot as plt
import seaborn

"采样点选取,0-2间取2^n个采样点,时域横坐标,单位s"
n=10
x=np.linspace(0,2,2**n)
# print(x)
"构造信号曲线,其中包含频率1Hz,2Hz,3Hz,4Hz"
y=2*np.sin(2*np.pi*1*x)/np.pi-2*np.sin(2*np.pi*2*x)/(2*np.pi)+2*np.sin(2*np.pi*3*x)/(3*np.pi)\\
  -2*np.sin(2*np.pi*4*x)/(4*np.pi)

plt.plot(x, y)
plt.show()
yy=fft(y) #快速傅里叶变换
yreal = yy.real  #获取实部
ymiag = yy.imag  #获取虚部
yf=abs(fft(y)) #取绝对值
yf1 = abs(fft(y))/len(x) #归一化处理
yf2 = yf1[range(int(len(x)/2))] ##对称性只取一半区间

"采样频率"
dt=x[1]-x[0]
fs = 1/dt
freqs = fs/2*np.linspace(0,1,int(2**n/2))
print(len(yf2))
print(len(freqs))
plt.subplot(211)
plt.plot(x, y)
plt.title(\'Original wave\')
plt.subplot(212)
plt.plot(freqs[0:10], yf2[0:10])
plt.title(\'FFT wave\')
plt.show()
快速傅里叶变换

从及结果可以看出,通过快速傅里叶变换,找到了四个阶次的频率,分别是1Hz,2Hz,3Hz,4Hz,与构造的信号函数频率一一对应。

 

  • 傅里叶变换处理发动机力(时域-频域)

    上述阐述已知信号源的傅里叶变换,其实信号源是否已知并不影响傅里叶变换的结果,构造已知信号的函数,仅仅是为了直观的理解这个变换输入输出。【这里并未对傅里叶变换的底层实现进行阐述,仅仅说明如何在发动机载荷分解过程中需要用到的功能】。

    回到上一章接种获取到的发动机曲轴中心位置的时间-力的数据。每一组时间-力(力矩)都如同上诉的信号源(只是未知信号)。

    

 

     

 

 

 

import numpy as np
from scipy.fftpack import  fft,ifft
import matplotlib.pyplot as plt
import seaborn
from scipy.interpolate import interp1d
"读取时间载荷"
data = np.loadtxt(r\'TimeLoad\\A750 _transpose.csv\',dtype=np.float, delimiter=",")
"获取时间采样点"
t=data[0]
"获取z向时间载荷"
z=data[3]
print(t)
print(z)
"原始数据采样,采样点数量按照2**n个进行,运用插值法,获取格则采样数据,"
n=13
NFFT = 2**n
x, delta_step = np.linspace(t[0], t[len(t) - 1], NFFT, endpoint=True, retstep=True)
y_curve = interp1d(t,z,kind=\'linear\')
y = y_curve(x)

plt.plot(x, y)
plt.show()
yy=fft(y) #快速傅里叶变换
yreal = yy.real  #获取实部
ymiag = yy.imag  #获取虚部
yf=abs(fft(y)) #取绝对值
yf1 = abs(fft(y))/len(x) #归一化处理
yf2 = yf1[range(int(len(x)/2))] ##对称性只取一半区间

"采样频率"
dt=x[1]-x[0]
fs = 1/dt
freqs = fs/2*np.linspace(0,1,int(2**n/2))
print(len(yf2))
print(len(freqs))
plt.subplot(211)
plt.plot(x, y)
plt.title(\'Original wave\')
plt.subplot(212)
plt.plot(freqs[0:100], yf2[0:100])
plt.title(\'FFT wave\')
plt.show()
z向载荷傅里叶变换

 

在怠速750r/min时,在25Hz的激励是62.157N。

 

以上是关于matlab中一数组的时域曲线如何转换为频域曲线的主要内容,如果未能解决你的问题,请参考以下文章

求助“利用MATLAB实现离散时间系统的时域和频域分析”的中译英翻译,数字信号处理,请高手帮帮忙

fft Matlab 振动时域到频率测量持续时间 = 60 s

如何对matlab plot生成的fig曲线图像进行去噪,平滑处理。

如何用matlab画平滑曲线?

第二章:整车发动机激励--快速傅里叶变换

曲线分类-特征提取