forked from dev-zzo/kiwiclient
-
Notifications
You must be signed in to change notification settings - Fork 35
/
proc_kiwi_iq_wav.m
120 lines (104 loc) · 3.44 KB
/
proc_kiwi_iq_wav.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
## Copyright (C) 2018 -*- octave -*-
## usage: [x,xx,fs,last_gpsfix]=proc_kiwi_iq_wav(wav_filename, max_last_gpsfix)
##
## input: wav_filename - filename of .wav file generated by kiwirecorder.py with options '-m iq --kiwi-wav'
## (optional) 2nd argument - maximal allowed age of the GNSS position solution (default is max_last_gpsfix=255)
## output: x - array of structs containing the information in the WAV chunks, i.e. IQ samples and GNSS timestamps
## xx - array of structs containing IQ samples and for each IQ samples a GNSS timestamp
## fs - sampling rate
## last_gpsfix - oldest GNSS position solution
##
## note: last_gps_fix = 255 means there was no GNSS position solution at all on the KiwiSDR
function [x,xx,fs,last_gpsfix]=proc_kiwi_iq_wav(fn, max_last_gpsfix)
if nargin == 0
print_usage;
end
if nargin == 1
max_last_gpsfix = 255;
end
[x,fs] = read_kiwi_iq_wav(fn);
if numel(x) < 2
xx = {};
last_gpsfix = 255;
return
endif
x = x(2:end); # omit the first chunk of 256 samples
last_gpsfix = max(cat(1,x.gpslast));
idx = find(cat(1,x.gpslast) < max_last_gpsfix);
## filter x
x = x(idx);
n = numel(x); ## from here on n is always >=2
## indices of fresh GNSS timestamps
idx_ts = find(diff(cat(1, x.gpslast)) < 0);
idx_ts = make_complete(idx_ts, n);
## remove outliers
ts = cat(1,x.gpssec) + 1e-9*cat(1,x.gpsnsec);
idx_ts = fix_outliers(ts, idx_ts, fs);
idx_ts = make_complete(idx_ts, n);
ts = ts(idx_ts);
## helper functions
gpssec = @(s) s.gpssec + 1e-9*s.gpsnsec;
get_fs = @(i) 512*(idx_ts(i) - idx_ts(i-1)) / (gpssec(x(idx_ts(i))) - gpssec(x(idx_ts(i-1))));
## keep history
fs = [];
## (1) interval before 1st fresh timestamp: extrapolate backwards
if idx_ts(1) > 1
fs(end+1) = get_fs(2);
k = 0;
for j=1:idx_ts(1)-1
xx(j).t = ts(1) - 512*(idx_ts(1)-1)/fs(end) + (512*k+[0:numel(x(j).z)-1])'/fs(end);
xx(j).z = x(j).z;
k += 1;
end
end
## (2) in between fresh timestamps: extrapolate between fresh timestamps
for i=2:numel(idx_ts)
##__fs = [ idx_ts(i) - idx_ts(i-1) gpssec(x(idx_ts(i))) - gpssec(x(idx_ts(i-1))) ]
fs(end+1) = get_fs(i);
k = 0;
for j=idx_ts(i-1):idx_ts(i)-1
xx(j).t = ts(i-1) + (512*k+[0:numel(x(j).z)-1])'/fs(end);
xx(j).z = x(j).z;
k += 1;
end
end
## (3) after the last fresh timestamp: extrapolate after the last fresh timestamp
k = 0;
for j=idx_ts(end):n
xx(j).t = ts(end) + (512*k+[0:numel(x(j).z)-1])'/fs(end);
xx(j).z = x(j).z;
k += 1;
end
if numel(fs) < 3 ## something wrong happened above and the sample rate value is not ok
last_gpsfix = 253;
end
## reduce the influence of outliers
fs = median(fs);
endfunction
function idx_ts=fix_outliers(t, idx_ts, fs)
## more than one sample difference -> outlier
outliers = find(abs(diff(t(idx_ts))*fs - 512*diff(idx_ts)) > 1);
b = ones(size(idx_ts))==1;
m = numel(outliers);
for i=1:2:m
if i+1 <= m
b(outliers(i)+1:outliers(i+1)+1) = false;
else
b(outliers(i)+1:end) = false;
end
end
## filter idx_ts
idx_ts(~b) = [];
endfunction
function idx_ts=make_complete(idx_ts, n)
## add fake fresh GNSS timestamps indices if numel(idx_ts)<2
if isempty(idx_ts)
idx_ts = [1 n];
elseif numel(idx_ts) == 1
if n-idx_ts > idx_ts
idx_ts = [idx_ts n];
else
idx_ts = [1 idx_ts];
end
end
endfunction