-
Notifications
You must be signed in to change notification settings - Fork 54
/
findPattern.m
executable file
·149 lines (142 loc) · 4.18 KB
/
findPattern.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
function [data, snr] = findPattern(N,OSR,NTF,ftest,Atest,f0,nlev,quadrature,debug)
% [data snr] = findPattern(N=1024,OSR=64,NTF,ftest,Atest,f0,nlev,quadrature,debug)
% findPattern.m Create a length-N bit-stream
% possessing good spectral properties when repeated.
%
% Iniitially created March 16 2000 for Jon Baldwin of ADI/CTS
% March 17 2000: Converted to function form, but not tested
% March 27 2001: Corrected some errors, generalized to bandpass, added two-tone support
% June 13 2003: Generalized to quadrature
%
%
% The procedure is:
% Try a bunch of signals. For each signal, calculate the noise power
% induced by resetting the state back to the value N samples earlier.
% This calculation is facilitated by pre-computation of the in-band
% portion of the FFTs of the impulse responses from each integrator's
% input to the output.
% Handle the input arguments
parameters = {'N' 'OSR' 'NTF' 'ftest' 'Atest' 'f0' 'nlev' 'quadrature' 'debug'};
defaults = { 2^13 64 NaN NaN 0.5 0 2 0 0 };
for i=1:length(defaults)
parameter = char(parameters(i));
if i>nargin | ( eval(['isnumeric(' parameter ') ']) & ...
eval(['any(isnan(' parameter ')) | isempty(' parameter ') ']) )
eval([parameter '=defaults{i};'])
end
end
if isnumeric(NTF) & isnan(NTF)
NTF = synthesizeNTF(6,OSR,1,(nlev+1)/2,f0);
end
if isnan(ftest)
if quadrature
if f0==0
ftest = 0.75/OSR; % Default ftest is 75% of pb width
else
ftest = f0 + 0.75/(2*OSR);
end
else
if f0==0
ftest = 0.75/(2*OSR); % Default ftest is 75% of pb width
else
ftest = f0 + 0.75/(4*OSR);
end
end
end
fbin = round(N*ftest);
ftest = fbin/N;
wtest = 2*pi*ftest;
Nsim = N + 2000; % +1000 is arbitrary
% Calculate the in-band portion of the FFTs of the impulse responses
if quadrature
ABCD = realizeQNTF(NTF,'FB');
[A,B,C,D] = partitionABCD(ABCD);
B2 = B(:,2);
order = size(A,1);
D = zeros(1,order);
A = A + B2*C;
B = eye(size(A));
sys = ss(mapQtoR(A), mapQtoR(B), mapQtoR(C), mapQtoR(D), 1);
imp = impulse(sys,N-1);
imp = squeeze( imp(:,1,1:2:2*order) +1i*imp(:,2,1:2:2*order) );
ABCD = mapQtoR(ABCD);
else
form = 'CRFB';
[a,g,b,c] = realizeNTF(NTF,form);
ABCD = stuffABCD(a,g,b,c,form);
[A,B,C,D] = partitionABCD(ABCD);
B2 = B(:,2);
order = size(A,1);
D = zeros(1,order);
A = A + B2*C;
B = eye(size(A));
sys = ss(A,B,C,D,1);
imp = impulse(sys,N-1);
imp = squeeze(imp);
end
if quadrature
if f0==0
ibb = [1:ceil(N/OSR+1)]; % in-band bins
else
ibb = [floor(N*(f0-1/(2*OSR))):ceil(N*(f0+1/(2*OSR)))] +1;
end
else
if f0==0
ibb = [1:ceil(N/(2*OSR)+1)]; % in-band bins
else
ibb = [max(1,floor(N*(f0-1/(4*OSR)))):ceil(N*(f0+1/(4*OSR)))] +1;
end
end
in_band = zeros(length(ibb),order);
for i = 1:order
spec = fft(imp(:,i))/(N/2);
in_band(:,i) = spec(ibb);
end
% Try a bunch of phases and slightly different amplitudes
Atest0 = Atest;
best_fom = Inf;
for i = 1:100 % 50 is arbitrary
Atest = Atest0 + 0.01*randn(1,1);
phi = 2*pi*rand(size(wtest));
u = zeros(1,Nsim+1);
for k = 1:length(wtest);
if quadrature
u = u + exp( 1i*( wtest(k)*[0:Nsim] + phi(k) ) );
else
u = u + sin( wtest(k)*[0:Nsim] + phi(k) );
end
end
u = u*Atest*(nlev-1)/length(wtest);
if quadrature
% [v xn] = simulateQDSM(u,ABCD,nlev);
u = [real(u); imag(u)];
[v xn] = simulateDSM(u,ABCD,nlev*[1 1] );
v = v(1,:) + 1i*v(2,:);
xn = xn(1:2:end,:) + 1i*xn(2:2:end,:);
else
[v xn] = simulateDSM(u,ABCD,nlev);
end
delta = xn(:,1:Nsim-N) - xn(:,N+1:Nsim);
spec = in_band * delta;
[fom start] = min( dbp( sum(abs(spec).^2) ) );
if fom < best_fom
best_fom = fom;
best = [Atest*(nlev-1) phi start];
data = v(start+1:start+N);
if debug
fprintf(1,'i=%d, fom=%g\n', i, fom);
end
end
end
snr = dbv(best(1)) - best_fom;
if 1 | debug
V0 = fft(data.*ds_hann(N))/(N*(nlev-1)/4);
snr0 = calculateSNR(V0(ibb),fbin-ibb(1)+1);
V1 = fft(data)/(N*(nlev-1)/2);
snr1 = -dbp( sum(abs(V1(ibb).^2)) / sum(abs(V1(fbin+1)).^2) -1 );
f = ibb-1;
plot( f,dbv(V0(ibb)), f,dbv(V1(ibb)) )
msg = sprintf('SNR = %.1fdB (%.1fdB)', snr1, snr0);
text(f(round(end/2)),-20,msg,'hor','center');
end
return