-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathhistogram_freq_domain.m
92 lines (71 loc) · 1.78 KB
/
histogram_freq_domain.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
% where the work is done
%% load data
addpath("data\")
load("data\data.mat")
%% frame + dft + store
next = 1;
OVERLAP_RATIO = 0.5;
FRAME_LENGTH = floor(0.02*fs);
S_cmp = [];
frame_count = 0;
while next<length(clean_1)
%% frame
[sl, next] = frame([clean_1;], next, ...
"overlap_ratio", OVERLAP_RATIO, ...
"length", FRAME_LENGTH);
frame_count = frame_count+1;
%% dft
Sl = fft(sl);
%% store
S_cmp(:,frame_count) = Sl;
%% processing
% Sl = Yl; % no processing
%% idft
% sl = ifft(Sl);
%% overlap
% output = attach_frame(output, sl, "overlap_ratio", OVERLAP_RATIO);
end
%% see histogram
K = 40;
S_sample = S_cmp(K,:);
index = ~(abs(S_sample) <= 1e-3);
S_sample = S_sample(index); % filter out zero-value samples (speech absent)
S_mag = abs(S_sample);
S_pha = phase(S_sample);
S_real = real(S_sample);
S_imag = imag(S_sample);
figure()
subplot(2,2,1)
histogram(S_mag, 250,...
"DisplayStyle", 'bar') %
subplot(2,2,2)
histogram(S_pha)
subplot(2,2,3)
histogram(S_real)
subplot(2,2,4)
histogram(S_imag)
%% deprecated
rayleigh_pd = fitdist(S_mag','rayleigh');
figure()
x_values = linspace(0,5,200);
y_values = pdf(rayleigh_pd, x_values);
plot(x_values, y_values)
%% histfit
figure()
pd_name = 'rayleigh';
histfit(S_mag',250,pd_name)
pdf = fitdist(S_mag',pd_name);
% log:
% birnbaumsaunders is a nice fit
% generalized pareto (gp) is also ok
% lognormal is ok
%% logpdf maximum likelihood estimate
custlogpdf = @(data,beta,gamma_p,nu) log(gamma_p)+gamma_p*log(beta)-...
log(gamma(nu))+(gamma_p*nu-1)*log(data)-beta*data.^gamma_p;
phat = mle(S_mag','logpdf',custlogpdf,'start',[0.5 2 1]);
%% get values
y=[];
for x = linspace(0,12,200)
y = [y, custlogpdf(x,phat(1),phat(2),phat(2))];
end
plot(x,exp(y))