forked from vicente-medina/TurbulenceToolbox
-
Notifications
You must be signed in to change notification settings - Fork 0
/
tt_spectrum_vic.m
60 lines (43 loc) · 1.31 KB
/
tt_spectrum_vic.m
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
function res = tt_spectrum_vic(VEL_IN)
%Number data
num_data = size(VEL_IN.x,1);
%Sampling data
Fs = 25;
%Number of samples per averaging
nfft = min(512,num_data);
%Data window
window = hanning(nfft);
%Overlaping
noverlap = 50;
%Power spectral density plot
%psd(VEL_IN.x,nfft,Fs,window,noverlap);
%Power spectral density data
[Pxx,f] = psd(VEL_IN.x,nfft,Fs,window,noverlap);
%Fit Kolmogorov
g = fittype('E*f^n','independent','f','dependent','Pxx','coefficients',{'E','n'});
%Options
opts = fitoptions(g);
opts.StartPoint = [100, -1];
%Fitting
fitted_line = fit(f(5:end),Pxx(5:end),g,opts)
Coef = coeffvalues(fitted_line);
%Equation
Eq = strcat('m_{Theory} = -1.6666 m_{Data} = ',num2str(1/Coef(2)));
%Line
YY=(f(5:end).^Coef(2))*Coef(1);
%Theory
YY2=(f(5:end).^(-3/5))*Coef(1);
%Plot
fig = figure;
set(fig,'Color',[1 1 1]);
loglog(f(5:end),Pxx(5:end));
line(f(5:end),YY,'LineWidth',2,'Color','black');
line(f(5:end),YY2,'LineWidth',2,'Color','red','LineStyle','--');
grid on;
xlabel('Frequency (s-1)');
ylabel('Energy (J·s)');
title(Eq);
Legend('Data','Fitting','Theory');
%Return
res = 0;
end