-
Notifications
You must be signed in to change notification settings - Fork 0
/
analiziraj.m
169 lines (139 loc) · 4.18 KB
/
analiziraj.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
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
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
% Časovna in spektralna analiza zvočnih posnetkov
function analiziraj(file1=0, file2=0)
pkg load signal;
global sf;
if (file1 != 0 && file2 != 0)
[x, fs] = audioread(file1);
[y, ~] = audioread(file2);
channels = audioinfo(file1).NumChannels;
else
global x y fs channels;
endif
save_graph = isfield(sf.param, "save_graphs_to_files") && ...
sf.param.save_graphs_to_files == true;
close all;
disp("-- Analiza signalov ---------");
tic;
% Če ima posnetek dva kanala, za izdelavo analize izračunaj njuno povprečje
if (channels == 2)
x_averaged = (x(:,1) + x(:,2))/2;
y_averaged = (y(:,1) + y(:,2))/2;
else
x_averaged = x;
y_averaged = y;
endif
if (numel(x)/fs > 60)
disp("# Omejujem dolžino posnetkov na 60 s pri FFT analizah");
x_reduced = x_averaged(1:60 * fs);
y_reduced = y_averaged(1:60 * fs);
else
x_reduced = x_averaged;
y_reduced = y_averaged;
endif
if (save_graph)
H1 = figure("visible", "off");
else
H1 = figure;
endif
pos = get(H1, 'Position');
set(H1, 'Position', [pos(1), pos(2), pos(3)*1.5, pos(4)]);
draw_amplitude(x_averaged, y_averaged, fs);
if (save_graph)
H2 = figure("visible", "off");
else
H2 = figure;
pos = get(H2, 'Position');
set(H2, 'Position', [pos(1), pos(2), pos(3)*1.5, pos(4)*1.5]);
subplot(2, 2, 1);
endif
draw_whole_spectrum(x_reduced, fs, "Močnostni spekter izvirnega posnetka");
if (save_graph)
H3 = figure("visible", "off");
else
subplot(2, 2, 2);
endif
draw_whole_spectrum(y_reduced, fs, "Močnostni spekter predelanega posnetka");
if (save_graph)
H4 = figure("visible", "off");
else
subplot(2, 2, 3);
endif
draw_spectrogram(x_reduced, fs, "Spektrogram izvirnega posnetka");
if (save_graph)
H5 = figure("visible", "off");
else
subplot(2, 2, 4);
endif
draw_spectrogram(y_reduced, fs, "Spektrogram predelanega posnetka");
printf("# čas analize = %.3f s\n", toc);
if (save_graph)
img_settings = ["-r80"];
disp("# shranjujem v slikovne datoteke");
print(H1, "-r90", "sf_ampl.png");
print(H2, img_settings, "sf_spekter_x.png");
print(H3, img_settings, "sf_spekter_y.png");
print(H4, img_settings, "sf_spektrogram_x.png");
print(H5, img_settings, "sf_spektrogram_y.png");
endif
endfunction
% Izris amplitude dveh posnetkov na en graf v odvisnosti od časa
function draw_amplitude(x, y, fs)
L = length(x);
t = (0:L-1) / fs;
plot(t, y);
hold on;
plot(t, x);
xlim([ t(1) t(end) ]);
ylim([ -1.5 1.5 ]);
ylabel("Rel. amplituda");
xlabel("Čas (s)");
title("Časovni potek amplitude");
grid on;
# set(gca, "fontsize", 12);
[~, hobj, ~, ~] = legend("Predelan", "Izvirni");
hl = findobj(hobj, "type", "line");
set(hl, "LineWidth", 2); % Odebeli črte v legendi
endfunction
% Fourierova analiza vseh vzorcev signala
function [ampl_max, freq_max] = draw_whole_spectrum(x, fs, name)
% Dolžini vzorcev
N = length(x);
Nfft = 2^nextpow2(N);
% Frekvenčni razpon z enosmerno komponento (0 Hz) v sredini
frequencies = (-Nfft/2:Nfft/2-1) * fs/Nfft;
% DFT signala (z oknjenjem?)
% Zamaknjen, da so negativne frek. na levi in pozitivne na desni
# amplitudes = abs(fftshift(fft(x, Nfft)) / N);
amplitudes = abs(fftshift(fft(x .* hamming(N), Nfft)) / N);
semilogy(frequencies/1000, amplitudes);
# plot(frequencies, 10 * log10(amplitudes));
xlim([ frequencies(1) frequencies(end) ]/1000);
xlabel("Frekvenca (kHz)");
ylabel("Amplituda (dB)");
title(name); % TODO annotation("textbox") ?
grid on;
[ampl_max, ampl_i] = max(10 * log10(amplitudes));
freq_max = frequencies(ampl_i);
# printf("%s: maksimum %.2f dB pri %.2f Hz\n", ...
# name, ampl_max, freq_max);
endfunction
% Kratkotrajna Fourierova analiza kosov signala (za izris spektrograma)
function draw_spectrogram(x, fs, name)
# N = 64; % Št. posameznih analiziranih kosov
step = ceil(50 * fs/1000);
window = ceil(100 * fs/1000);
N = 2^nextpow2(window);
[S, f, t] = specgram(x, N, fs, window, window - step);
imagesc(t, f/1000, 20*log10(abs(S)));
set(gca(), "ydir", "normal");
xlabel("Čas (s)");
ylabel("Frekvenca (kHz)");
title(name);
endfunction
# try
# clf(2);
# catch
# pos = get(gcf,'Position');
# % x, y, širina, višina
# set(gcf, 'Position', [pos(1), pos(2), pos(3)*2, pos(4)*2]);
# end_try_catch