Matlab - PSD de una función temporal

 
Vista:
sin imagen de perfil
Val: 3
Ha aumentado su posición en 24 puestos en Matlab (en relación al último mes)
Gráfica de Matlab

PSD de una función temporal

Publicado por Daniel (2 intervenciones) el 27/04/2019 19:15:30
Hola a todos,

Dispongo de los datos de aceleración recogidos de un acelerómetro, y del tiempo de adquisición, siendo la gráfica aceleración vs tiempo la que muestro en la siguiente imagen:

1

Mi duda es, a partir de esta gráfica, cómo podría calcular en Matlab la función PSD en el dominio de la frecuencia, a partir de la función temporal arriba mostrada?

A modo de ejemplo del problema práctico, he empleado esta función para generar la gráfica mostrada:

1
2
3
t = 0:.001:0.25;
x = sin(2*pi*50*t) + sin(2*pi*120*t);
plot(t,x)


Muchas gracias por la ayuda!
Valora esta pregunta
Me gusta: Está pregunta es útil y esta claraNo me gusta: Está pregunta no esta clara o no es útil
0
Responder
Imágen de perfil de David Correa
Val: 497
Ha mantenido su posición en Matlab (en relación al último mes)
Gráfica de Matlab

PSD de una función temporal

Publicado por David Correa (1094 intervenciones) el 28/04/2019 12:44:34
Hola Daniel;

Muy interesante lo que deseas hacer.

En algún momento trabaje con datos de acelerometros para evaluar variaciones de movimiento y evalue varios aspectos que podrían ser de tu interes, a continuación un ejemplo:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
clear all, close all,clc, warning off
 
fname = ['accelerometer_copia.xls'];
[dat,txt,all] = xlsread(fname);
 
t = dat(:,1);
x = dat(:,2);
y = dat(:,3);
z = dat(:,4);
 
nt = [100:100:max(t)];
nx = interp1(t,x,nt);
ny = interp1(t,y,nt);
nz = interp1(t,z,nt);
 
clear dat all txt fname
 
%%
figure(1)
subplot(3,1,1)
plot((nt./1000)./60,nx)
ylabel('X [m/s2]')
xlabel('Tiempo [min]')
grid on
 
subplot(3,1,2)
plot((nt./1000)./60,ny)
ylabel('Y [m/s2]')
xlabel('Tiempo [min]')
grid on
 
subplot(3,1,3)
plot((nt./1000)./60,nz)
ylabel('Z [m/s2]')
xlabel('Tiempo [min]')
grid on
 
 
%%
Fs = 1000*(1./100);      % 10 oscilaciones x segundo
 
L = length(nt);
NFFT = 4096;
X = fft(nx,NFFT)/L;     X(1:2) = nan;
Y = fft(ny,NFFT)/L;     Y(1:2) = nan;
Z = fft(nz,NFFT)/L;     Z(1:2) = nan;
 
f = Fs/2*linspace(0,1,NFFT/2+1);
 
% Plot de un lado del espectro de amplitud.
figure
subplot(311)
plot(f,2*abs(X(1:NFFT/2+1)))
title('Espectro de Amplitud X(t)')
xlabel('Frecuencia (Hz)')
ylabel('|X(f)|')
% ylim([0 0.04])
grid on
 
 
subplot(312)
plot(f,2*abs(Y(1:NFFT/2+1)))
title('Espectro de Amplitud Y(t)')
xlabel('Frecuencia (Hz)')
ylabel('|Y(f)|')
% ylim([0 0.04])
grid on
 
 
subplot(313)
plot(f,2*abs(Z(1:NFFT/2+1)))
title('Espectro de Amplitud Z(t)')
xlabel('Frecuencia (Hz)')
ylabel('|Z(f)|')
% ylim([0 0.04])
grid on
 
 
 
%%
figure
subplot(311)
plot(1./f,2*abs(X(1:NFFT/2+1)))
title('Espectro de Amplitud X(t)')
xlabel('Periodo (Seg)')
ylabel('|X(f)|')
% ylim([0 0.04])
grid on
 
 
subplot(312)
plot(1./f,2*abs(Y(1:NFFT/2+1)))
title('Espectro de Amplitud Y(t)')
xlabel('Periodo (Seg)')
ylabel('|Y(f)|')
% ylim([0 0.04])
grid on
 
 
subplot(313)
plot(1./f,2*abs(Z(1:NFFT/2+1)))
title('Espectro de Amplitud Z(t)')
xlabel('Periodo (Seg)')
ylabel('|Z(f)|')
% ylim([0 0.04])
grid on
 
 
 
%%
 
figure
subplot(311)
[a,b] = hist(nx,[-1:0.1:12]);
bar(b,100*a/sum(a),1)
xlabel('Acelerometro X [m/s2]')
grid on
set(gca,'xtick',[-1:0.5:12])
xlim([-1 12])
ylabel('%')
 
subplot(312)
[a,b] = hist(ny,[-1:0.1:12]);
bar(b,100*a/sum(a),1)
xlabel('Acelerometro Y [m/s2]')
grid on
set(gca,'xtick',[-1:0.5:12])
xlim([-1 12])
ylabel('%')
 
subplot(313)
[a,b] = hist(nz,[-1:0.1:12]);
bar(b,100*a/sum(a),1)
xlabel('Acelerometro Z [m/s2]')
grid on
set(gca,'xtick',[-1:0.5:12])
xlim([-1 12])
ylabel('%')

Te adjunto los archivos.

​Saludos
David Correa
Director de Servicios de Programación
E-mail: [email protected]
Web page: https://www.fismatlab.com
Facebook: https://www.facebook.com/fismatlabperu
Blog: http://fismatlab.blogspot.com
Spot: https://www.youtube.com/watch?v=NTDY-MRnFMk
WhatsApp: +51 - 922210488
Valora esta respuesta
Me gusta: Está respuesta es útil y esta claraNo me gusta: Está respuesta no esta clara o no es útil
0
Comentar
sin imagen de perfil
Val: 3
Ha aumentado su posición en 24 puestos en Matlab (en relación al último mes)
Gráfica de Matlab

PSD de una función temporal

Publicado por Daniel (2 intervenciones) el 29/04/2019 14:09:35
Muchas gracias por tu respuesta, muy interesante y me ha aclarado varios conceptos.

Mi post viene dado ya que estoy realizando un análisis de vibración aleatoria ("Random Vibration" en Ansys Workbench) de un componente de automoción, el cuál ha sido monitorizado con acelerómetros durante una ruta de carretera para obtener las cargas.

Dicho esto, en el análisis Random Vibration, las cargas hay que introducirlas como una densidad de potencia espectral (PSD), que vendría a ser una representación, en este caso, de la aceleración al cuadrado por Hz, frente a frecuencia, es decir: ((m/s2)^2)/Hz - Hz. Un ejemplo podría ser el que muestro en la siguiente imagen:


2

Mi duda sería pues, cómo calcular la gráfica PSD mostrada en la imagen derecha.

Muchas gracias.
Valora esta respuesta
Me gusta: Está respuesta es útil y esta claraNo me gusta: Está respuesta no esta clara o no es útil
0
Comentar
Imágen de perfil de David Correa
Val: 497
Ha mantenido su posición en Matlab (en relación al último mes)
Gráfica de Matlab

PSD de una función temporal

Publicado por David Correa (1094 intervenciones) el 30/04/2019 13:07:34
Hola Daniel;

Te sugiero que utilices la transformada de Fourier, a través de esa función podrá hallar el espectro de potencia.

​Saludos
David Correa
Director de Servicios de Programación
E-mail: [email protected]
Web page: https://www.fismatlab.com
Facebook: https://www.facebook.com/fismatlabperu
Blog: http://fismatlab.blogspot.com
Spot: https://www.youtube.com/watch?v=NTDY-MRnFMk
WhatsApp: +51 - 922210488
Valora esta respuesta
Me gusta: Está respuesta es útil y esta claraNo me gusta: Está respuesta no esta clara o no es útil
0
Comentar