summaryrefslogtreecommitdiffstats
path: root/SD-VBS/benchmarks/texture_synthesis/src/matlab/modkurt.m
diff options
context:
space:
mode:
Diffstat (limited to 'SD-VBS/benchmarks/texture_synthesis/src/matlab/modkurt.m')
-rwxr-xr-xSD-VBS/benchmarks/texture_synthesis/src/matlab/modkurt.m172
1 files changed, 0 insertions, 172 deletions
diff --git a/SD-VBS/benchmarks/texture_synthesis/src/matlab/modkurt.m b/SD-VBS/benchmarks/texture_synthesis/src/matlab/modkurt.m
deleted file mode 100755
index e144d32..0000000
--- a/SD-VBS/benchmarks/texture_synthesis/src/matlab/modkurt.m
+++ /dev/null
@@ -1,172 +0,0 @@
1function [chm, snrk] = modkurt(ch,k,p);
2
3% Modify the kurtosis in one step, by moving in gradient direction until
4% reaching the desired kurtosis value.
5% It does not affect the mean nor the variance, but it affects the skewness.
6% This operation is not an orthogonal projection, but the projection angle is
7% near pi/2 when k is close to the original kurtosis, which is a realistic assumption
8% when doing iterative projections in a pyramid, for example (small corrections
9% to the channels' statistics).
10%
11% [chm, snrk] = modkurt(ch,k,p);
12% ch: channel
13% k: desired kurtosis (k=M4/M2^2)
14% p [OPTIONAL]: mixing proportion between k0 and k
15% it imposes (1-p)*k0 + p*k,
16% being k0 the current kurtosis.
17% DEFAULT: p = 1;
18
19
20% Javier Portilla, Oct.12/97, NYU
21
22Warn = 0; % Set to 1 if you want to see warning messages
23if ~exist('p'),
24 p = 1;
25end
26
27me=mean2(ch);
28ch=ch-me;
29
30% Compute the moments
31
32m=zeros(12,1);
33for n=2:12,
34 m(n)=mean2(ch.^n);
35end
36
37% The original kurtosis
38
39k0=m(4)/m(2)^2;
40snrk = snr(k, k-k0);
41if snrk > 60,
42 chm = ch+me;
43 return
44end
45k = k0*(1-p) + k*p;
46
47% Some auxiliar variables
48
49a=m(4)/m(2);
50
51% Coeficients of the numerator (A*lam^4+B*lam^3+C*lam^2+D*lam+E)
52
53A=m(12)-4*a*m(10)-4*m(3)*m(9)+6*a^2*m(8)+12*a*m(3)*m(7)+6*m(3)^2*m(6)-...
54 4*a^3*m(6)-12*a^2*m(3)*m(5)+a^4*m(4)-12*a*m(3)^2*m(4)+...
55 4*a^3*m(3)^2+6*a^2*m(3)^2*m(2)-3*m(3)^4;
56B=4*(m(10)-3*a*m(8)-3*m(3)*m(7)+3*a^2*m(6)+6*a*m(3)*m(5)+3*m(3)^2*m(4)-...
57 a^3*m(4)-3*a^2*m(3)^2-3*m(4)*m(3)^2);
58C=6*(m(8)-2*a*m(6)-2*m(3)*m(5)+a^2*m(4)+2*a*m(3)^2+m(3)^2*m(2));
59D=4*(m(6)-a^2*m(2)-m(3)^2);
60E=m(4);
61
62% Define the coefficients of the denominator (F*lam^2+G)^2
63
64F=D/4;
65G=m(2);
66
67% test
68test = 0;
69
70if test,
71
72 grd = ch.^3 - a*ch - m(3);
73 lam = -0.001:0.00001:0.001;
74 k = (A*lam.^4+B*lam.^3+C*lam.^2+D*lam+E)./...
75 (F*lam.^2 + G).^2;
76 for lam = -0.001:0.00001:0.001,
77 n = lam*100000+101;
78 chp = ch + lam*grd;
79 k2(n) = mean2(chp.^4)/mean2(chp.^2)^2;
80 %k2(n) = mean2(chp.^4);
81 end
82 lam = -0.001:0.00001:0.001;
83 snr(k2, k-k2)
84
85end % test
86
87% Now I compute its derivative with respect to lambda
88% (only the roots of derivative = 0 )
89
90d(1) = B*F;
91d(2) = 2*C*F - 4*A*G;
92d(3) = 4*F*D -3*B*G - D*F;
93d(4) = 4*F*E - 2*C*G;
94d(5) = -D*G;
95
96mMlambda = roots(d);
97
98tg = imag(mMlambda)./real(mMlambda);
99mMlambda = mMlambda(find(abs(tg)<1e-6));
100lNeg = mMlambda(find(mMlambda<0));
101if length(lNeg)==0,
102 lNeg = -1/eps;
103end
104lPos = mMlambda(find(mMlambda>=0));
105if length(lPos)==0,
106 lPos = 1/eps;
107end
108lmi = max(lNeg);
109lma = min(lPos);
110
111lam = [lmi lma];
112mMnewKt = polyval([A B C D E],lam)./(polyval([F 0 G],lam)).^2;
113kmin = min(mMnewKt);
114kmax = max(mMnewKt);
115
116% Given a desired kurtosis, solves for lambda
117
118if k<=kmin & Warn,
119 lam = lmi;
120 warning('Saturating (down) kurtosis!');
121 kmin
122elseif k>=kmax & Warn,
123 lam = lma;
124 warning('Saturating (up) kurtosis!');
125 kmax
126else
127
128% Coeficients of the algebraic equation
129
130c0 = E - k*G^2;
131c1 = D;
132c2 = C - 2*k*F*G;
133c3 = B;
134c4 = A - k*F^2;
135
136% Solves the equation
137
138r=roots([c4 c3 c2 c1 c0]);
139
140% Chose the real solution with minimum absolute value with the rigth sign
141
142tg = imag(r)./real(r);
143%lambda = real(r(find(abs(tg)<1e-6)));
144lambda = real(r(find(abs(tg)==0)));
145if length(lambda)>0,
146 lam = lambda(find(abs(lambda)==min(abs(lambda))));
147 lam = lam(1);
148else
149 lam = 0;
150end
151
152end % if ... else
153
154
155% Modify the channel
156
157chm=ch+lam*(ch.^3-a*ch-m(3)); % adjust the kurtosis
158chm=chm*sqrt(m(2)/mean2(chm.^2)); % adjust the variance
159chm=chm+me; % adjust the mean
160
161% Check the result
162%k2=mean2((chm-me).^4)/(mean2((chm-me).^2))^2;
163%SNR=snr(k,k-k2)
164
165
166
167
168
169
170
171
172