BP.zip_BP预测_信号预测

% bpfig34: BP Figure 34 -- Phase plane evolution at % BP-Interior iterations. % % BP-Interior iteratively sparsifies the MOF solution. % % Signal: FM-Cosine % Signal Length: 1024 % Dictionary: Cosine Packets with D = log2(n)-4 % Observations: % (a) Initial phase plane, i.e. the MOF phase plane. % (b) Phase plane at the 1st BP-Interior iteration. % (c) Phase plane at the 2st BP-Interior iteration. % (d) Phase plane at the 3st BP-Interior iteration. % (e) Phase plane at the 4st BP-Interior iteration. % (f) Phase plane at the BP-Interior iteration teramination. % help bpfig34 n = 1024; D = log2(n); qmf = MakeONFilter('Symmlet', 8);bell = 'Sine'; x = InputSignal('FM', n); cBPMovie= zeros(n * (D-4+1), 6); NameOfDict = 'CP'; par1 = log2(n) - 4; par2 = 0; par3 = 0; MaxBPIter = 30; CGAccuracy = 1e-1; FeaTol = 1e-1; OptTol = 1e-1; %*** Size of the problem b = x(:); n = length(b); [m L] = SizeOfDict(n, NameOfDict, par1, par2, par3); %*** Setting the controlling parameters delta = OptTol; gamma = OptTol; sigma = .995; MaxCGIter = m * log2(m); %*** Initializing (x, y, z), make sure x and z are INTERIOR %*** Initializing by MOF representation fprintf('\nBP-Interior:\n\n'); disp('Initializing BP-Interior by MOF ...') x = MOF(b, NameOfDict, par1, par2, par3,1e-5); x = [ x .* (x > 0) ; (-x) .* (x < 0) ]; sg = sign(x); y = FastSS(sg, n, NameOfDict, par1, par2, par3); ap = FastAA(y,NameOfDict,par1,par2, par3);ap = ap(1:m); mp = 1.1 * max(max(abs(ap))); ap = ap ./ mp; y = y ./ mp; z1 = 1- ap(:); z2 = 1+ ap(:); z = [z1 ; z2]; x = x + .1; x0 = x; y0 = y; z0 = z; %*** Initialization for BP-Interior Iterations mu0 = .01; f = [z .* x; b - FastSS(x, n, NameOfDict, par1, par2, par3) - delta * y; ... - FastAA(y, NameOfDict, par1, par2, par3) - z + gamma * x + 1]; mu = mu0 * norm(f) / sqrt(2 * m); k = -1; gamma = gamma ^ 2; delta = delta ^ 2; zerosn = zeros(n,1); %*** Iteration disp(' ') disp('BP-Interior Iterations ...') head = ' mItn PrimalFea DualFea DualGap Obj CGStep'; disp(head) for i = 1:MaxBPIter, if i >= 1 & i <= 5, cBPMovie(:,i) = x(1:m)-x((m+1):(2*m)); end % Calculate the Residules v1 = mu - z .* x; v2 = b - FastSS(x, n, NameOfDict, par1, par2, par3); v3 = - FastAA(y, NameOfDict, par1, par2, par3) - z + 1; % Test Termination PrimalFeas = norm(v2) / (1 + norm(x)); DualFeas = norm(v3) / (1 + norm(y)); PDGap = z' * x / (1 + abs(sum(x))); Obj = sum(abs(x)); fprintf('BPStep %4g: %8.2e %8.2e %8.2e %14.7e %3g\n', i-1, PrimalFeas, DualFeas, PDGap, Obj, k) Feasible = PrimalFeas < FeaTol & DualFeas < FeaTol; Converged = PDGap < OptTol; if Feasible & Converged, break; end % Solve Newton's Direction by CG: (SHA + delta I) * dy = h H = 1 ./ (z ./ x + gamma); h = v2 - FastSS(H .* (v1 ./ x - v3), n, NameOfDict, par1, par2, par3); thresh = CGAccuracy/i * norm(h); thresh = thresh ^ 2; dy = zerosn; r = h; rho_1 = norm(r) ^ 2; for k = 1:MaxCGIter, %fprintf(' CG Step: %3g: %g\n', k, rho_1); if k == 1, p = r; else beta = rho_1 / rho_0; p = r + beta * p; end w = FastSS(H .* FastAA(p, NameOfDict, par1, par2, par3), n, ... NameOfDict, par1, par2, par3) + delta * p; alpha = rho_1 / (p' * w); r = r - alpha * w; rho_0 = rho_1; rho_1 = norm(r) ^2; dy = dy + alpha * p; if rho_1 < thresh, break; end end if k == MaxCGIter, if imp < 0, x = xold; end fprintf('Too many CG iterations!\n') break; end dx=H .* (v1 ./ x - v3) + H .* FastAA(dy, NameOfDict, par1, par2, par3); dz = v1 ./ x - dx .* z ./ x; negdx = dx < -1e-14; negdz = dz < -1e-14; Rp = min(x(negdx) ./ abs(dx(negdx))); Rd = min(z(negdz) ./ abs(dz(negdz))); % Update (x y z mu) xold = x; x = x + sigma * Rp * dx; y = y + sigma * Rd * dy; z = z + sigma * Rd * dz; mu = (1 - min([Rp Rd sigma])) * mu; end %% Feasibility: %% Let dx = x - x0 %% Project dx onto the null space of \Phi by Conjugate Gradients: %% dx = (I - \Phi^T * (\Phi * \Phi^T)^-1 * \Phi) dx %% x = x0 + dx %% dx = x - x0; r = FastSS(dx, n, NameOfDict, par1, par2, par3); thresh = 1e-16 * norm(r); thresh = thresh ^ 2; dd = zerosn; rho_1 = norm(r) ^ 2; for k = 1:MaxCGIter, %fprintf(' CG Step: %3g: %g\n', k, rho_1); if k == 1, p = r; else beta = rho_1 / rho_0; p = r + beta * p; end w = FastSS(FastAA(p, NameOfDict, par1, par2, par3), n, ... NameOfDict, par1, par2, par3); alpha = rho_1 / (p' * w); r = r - alpha * w; rho_0 = rho_1; rho_1 = norm(r) ^2; dd = dd + alpha * p; if rho_1 < thresh, break; end end dx = dx - FastAA(dd, NameOfDict, par1, par2, par3); %dx = dx - FastAA(FastSS(dx, n, NameOfDict, par1, par2, par3), ... % NameOfDict, par1, par2, par3) / (2 * L); x = x0 + dx; c = x(1:m) - x((m+1):(2*m)); cBPMovie(:,6) = c; fprintf('Obj = %14.7e\n', sum(abs(c))); TFScale = 20; figure(1);clg for i = 1:6, subplot(3,2,i); PhasePlane(cBPMovie(:, i), 'CP', n, 0, 256, TFScale); end subplot(3,2,1);title('Phase Plane: BP Iteration = 0'); subplot(3,2,2);title('Phase Plane: BP Iteration = 1'); subplot(3,2,3);title('Phase Plane: BP Iteration = 2'); subplot(3,2,4);title('Phase Plane: BP Iteration = 3'); subplot(3,2,5);title('Phase Plane: BP Iteration = 4'); subplot(3,2,6);title('Phase Plane: BP Termination');