diff --git a/common/alpha.m b/common/alpha.m new file mode 100755 index 0000000..8936be9 --- /dev/null +++ b/common/alpha.m @@ -0,0 +1,13 @@ +% a = alpha(xs,xl,nl,pl,es,dts,pa) +function a = alpha(xs,xl,nl,pl,es,dts,pa) + +N = lge(xs); + +for n = 1:N, + if((xs(n) > pl*xl(n)) & (dts(n)*es(n) >= pa*es(n))) +% a(n) = ((es(n)-nl(n)+0.001)/(es(n)+0.001))^2; + a(n) = max(1-nl(n)/(es(n)+0.008),0)^2; + else + a(n) = 0; + end; +end; diff --git a/common/ccorr.m b/common/ccorr.m new file mode 100755 index 0000000..019932a --- /dev/null +++ b/common/ccorr.m @@ -0,0 +1,52 @@ +% Kreuzkorrelation von x und y +% ---------------------------- +% +% Aufruf: +% Function [rxy,fpo]=ccorr(x,y,N) +% +% Parameter: +% x : Eingangsdaten 1 +% y : Eingangsdaten 2 +% N : Maximale Verschiebung der Eingangsdaten +% +% Rückgabewerte: +% rxy: Kreuzkorrelationskoeffzienten +% fpo: Anzahl der Floating-Point Operationen + +% ------------------------------------------------------------------------- +% Datum : 07.08.2001 +% Autor : Jens Ahrensfeld +% Thema : Diplomarbeit +% Datei : ccorr.m +% Benötigte Dateien: +% +% ------------------------------------------------------------------------- +function [rxy,fpo]=ccorr(x,y,N) + +% Ermittling der FFT-Länge (für Radix-2) +C = 2^ceil(log2(N)); + +% Zero-Padding wenn N != 2^r +xzp = [x;zeros(C-N,1)]; +yzp = [y;zeros(C-N,1)]; + +flops(0); + +% C-Punkte FFT der Eingangsdaten +Xk = fft(xzp,C); +Yk = fft(yzp,C); + +% Korrelation im Frequenzbereich +Rk = Xk.*conj(Yk); + +% Rücktransformation +rk = ifft(Rk,C)/C; + +% Abspeichern der ersten N Werte +rxy(1:N/2) = rk(N/2+1:N); +rxy(N/2+1:N) = rk(1:N/2); + +fpo = flops; + +% ------------------------------------------------------------------------- +% Ende ccorr.m diff --git a/common/ccorr2.m b/common/ccorr2.m new file mode 100755 index 0000000..e43e0e0 --- /dev/null +++ b/common/ccorr2.m @@ -0,0 +1,32 @@ +% Kreuzkorrelation von x und y +% ---------------------------- +% +% Aufruf: +% Function [rxy]=ccorr2(x,y,M) +% +% Parameter: +% +% Rückgabewerte: + +% ------------------------------------------------------------------------- +% Datum : ..2001 +% Autor : Jens Ahrensfeld +% Thema : Diplomarbeit +% Datei : newfunc.m +% Benötigte Dateien: +% +% ------------------------------------------------------------------------- +function [rxy]=ccorr2(x,y,M) + +% Berechne Erwartungswert xm, ym von X und Y im Intervall N +for n = N:N:L-W, + rmax(k) = 0; + for w = 0:W-1 + sxy(k) = sum(x(n+w:-1:n-N+1+w).*y(n:-1:n-N+1)); + sxym(k) = sum(abs(x(n+w:-1:n-N+1+w).*y(n:-1:n-N+1))); + rmax(k) = max(sxy(k)/sxym(k),rmax(k)); + end; +% r(n:-1:n-N+1) = ones(N,1)*sxy(k)/sxym(k); + r(n:-1:n-N+1) = ones(N,1)*min(rmax(k),1); + k = k+1; +end; diff --git a/common/cdkap2.m b/common/cdkap2.m new file mode 100755 index 0000000..5d6e822 --- /dev/null +++ b/common/cdkap2.m @@ -0,0 +1,2 @@ +cd f: +cd \work\dpa\matlab\sim\kap2 \ No newline at end of file diff --git a/common/cdkap4.m b/common/cdkap4.m new file mode 100755 index 0000000..0bee277 --- /dev/null +++ b/common/cdkap4.m @@ -0,0 +1,2 @@ +cd f: +cd \work\dpa\matlab\sim\kap4 \ No newline at end of file diff --git a/common/cic_filter.m b/common/cic_filter.m new file mode 100755 index 0000000..02666e1 --- /dev/null +++ b/common/cic_filter.m @@ -0,0 +1,30 @@ +function [y, gain] = cic_filter(R, M, N, x) +% [y, gain] = cic_filter(R, M, N, x) + +a_i = [1 -1]; +b_i = 1; + +a_c = 1; +b_c = [1 zeros(1, M-1) -1]; + +gain = (M*R)^N +k_i = 1/(M*R); + +i_out = x; + +% Integrator +for i=1:N, + i_out = k_i*filter(b_i, a_i, i_out); +end; + +% Decimate +c_out = i_out(R:R:length(i_out)); + +% Comb +for i=1:N, + c_out = filter(b_c, a_c, c_out); +end; + + +y = c_out; + diff --git a/common/cielab2xyz.m b/common/cielab2xyz.m new file mode 100755 index 0000000..c91b246 --- /dev/null +++ b/common/cielab2xyz.m @@ -0,0 +1,45 @@ +% function [xyz] = cielab2xyz(lab, ref_xyz) +% +% // C-Code +% var_Y = ( CIE-L* + 16 ) / 116 +% var_X = CIE-a* / 500 + var_Y +% var_Z = var_Y - CIE-b* / 200 +% +% if ( var_Y^3 > 0.008856 ) var_Y = var_Y^3 +% else var_Y = ( var_Y - 16 / 116 ) / 7.787 +% if ( var_X^3 > 0.008856 ) var_X = var_X^3 +% else var_X = ( var_X - 16 / 116 ) / 7.787 +% if ( var_Z^3 > 0.008856 ) var_Z = var_Z^3 +% else var_Z = ( var_Z - 16 / 116 ) / 7.787 +% +% X = ref_X * var_X //ref_X = 95.047 Observer= 2°, Illuminant= D65 +% Y = ref_Y * var_Y //ref_Y = 100.000 +% Z = ref_Z * var_Z //ref_Z = 108.883 +% +% Source http://www.easyrgb.com + +function [xyz] = cielab2xyz(lab, ref_xyz) + +[m,n] = size(lab); +if ((m ~= 3) & (n == 3)) + lab = lab' +end; + +xyz_k(2,:) = (lab(1,:) + 16)/116; +xyz_k(1,:) = lab(2,:)/500 + xyz_k(2,:); +xyz_k(3,:) = xyz_k(2,:) - lab(3,:)/200; + +xyz = f1(xyz_k); + +xyz(1,:) = xyz(1,:).*ref_xyz(1,:); +xyz(2,:) = xyz(2,:).*ref_xyz(2,:); +xyz(3,:) = xyz(3,:).*ref_xyz(3,:); + +function xyz = f1(xyz_k) +thresh_n = 0.008856; + +[idx_m, idx_n] = find(xyz_k.^3 > thresh_n); +xyz(idx_m,idx_n) = xyz_k(idx_m,idx_n).^3; + +[idx_m, idx_n] = find(xyz_k.^3 <= thresh_n); +xyz(idx_m,idx_n) = (xyz_k(idx_m,idx_n) - 16/116)/7.787; diff --git a/common/clusters.m b/common/clusters.m new file mode 100755 index 0000000..f864d57 --- /dev/null +++ b/common/clusters.m @@ -0,0 +1,29 @@ +% clusters.m +% function cb = clusters(input,ref,maxIt) +% clusters(randn(2,400)/sqrt(12),randn(2,1)/sqrt(12)); +% 17.01.2003, Jens Ahrensfeld + +function cb = clusters(input,ref,maxIt) +eps = 0.4; +epse = 0.001; + +[nDim, nData] = size(input); +[nDimr, nDatar] = size(ref); + +if (nDim ~= nDimr), + error('Dimensions must match!'); +end; + +d = zeros(nDatar,1); + +for i=1:maxIt, + n = round(rand(1,1)*(nData-1))+1; % Randomly pick input vector + for k=1:nDatar, + d(k) = sqrt(sum((input(:,n)-ref(:,k)).^2)); % Calc distances between references and input + end; + [min, winidx] = min(d); % Calc best reference vector - the winner + winner = ref(:,winidx); + epsi = eps*(1-i/maxIt)+epse; + ref(:,winidx) = winner + epsi*(input(:,n)-winner); % Move winner towards input +end; +cb = ref; diff --git a/common/clusters_ob.m b/common/clusters_ob.m new file mode 100755 index 0000000..cf1c701 --- /dev/null +++ b/common/clusters_ob.m @@ -0,0 +1,144 @@ +% cltest.m +% function cltest(dim, numvecperclass, classdist, numClasses, numref, maxIt); +% Sample cltest(2, 80, 8, 3, 3, 800) +% 17.01.2003, Jens Ahrensfeld +% From +% "Color Quantization of Images", Orchard, Bouman +% IEEE Trans. on Sig. Proc., vol. 39, no. 12, pp. 2677-2690, Dec. 1991. + +function [cb, idx] = clusters_ob(xs,M,tse_max) +[nc,len] = size(xs); +% Root node +tree(1) = node_root(xs, 1:len); + +tse = 0; +node = tree(1); +end_nodes = 1; +num_nodes = 1; +lambdas = tree(1).lambda; +map = node.q; +for m=1:M-1, + + if (num_nodes == 0) + break; + end; + + [v, i] = max(lambdas); + n = end_nodes(i); + + node = tree(n); + + L_n = node.lambda; + num_nodes = num_nodes - 1; + + if(node.tse < tse_max) + break; + end; + + if L_n == 0 + continue; + end; + + lambdas(i) = []; + end_nodes(i) = []; + + [node2n, node2n1] = node_split(xs, node); + + if((node2n.N * node2n1.N) == 0) + disp('No further split'); + continue; + end; + + tree(2*n) = node2n; + tree(2*n+1) = node2n1; + + L2n = tree(2*n).lambda; + L2n1 = tree(2*n+1).lambda; + + end_nodes = [end_nodes 2*n 2*n+1]; + lambdas = [lambdas L2n L2n1]; + + num_nodes = num_nodes + 2; + + tse = tse + tree(n).tse; + + +end; +map = [tree(end_nodes).q]; +lambda = [tree(end_nodes).lambda]; +len = length(lambda); +[v,i] = sort(lambda); +s = i(len:-1:1); +cnt = 1; +for s=s, + if (isempty(tree(s).C)) + break + end + idx{cnt} = tree(s).C; + cb(:,cnt) = map(:,s); + cnt = cnt + 1; +end; +v; + +% --------------------------------------------------------- +function [node_2n, node_2n1] = node_split(xs, node) +C_n = node.C; +e_n = node.e; +q_n = node.q; + +i = find(e_n'*xs(:,C_n) <= e_n'*q_n); +C_2n = C_n(i); +i = find(e_n'*xs(:,C_n) > e_n'*q_n); +C_2n1 = C_n(i); + +xs_2n = xs(:,C_2n); +xs_2n1 = xs(:,C_2n1); + +R_n = node.R; +m_n = node.m; +N_n = node.N; + +R_2n = xs_2n*xs_2n'; +m_2n = sum(xs_2n,2); +N_2n = length(C_2n); + +R_2n1 = R_n - R_2n; +m_2n1 = m_n - m_2n; +N_2n1 = N_n - N_2n; + +node_2n = struct('R',R_2n,'m',m_2n,'N',N_2n,'C',C_2n,'q',0,'e',0,'lambda',0, 'tse', 0); +node_2n1 = struct('R',R_2n1,'m',m_2n1,'N',N_2n1,'C',C_2n1,'q',0,'e',0,'lambda',0, 'tse', 0); + +node_2n = node_stat(node_2n, xs); +node_2n1 = node_stat(node_2n1, xs); + +% --------------------------------------------------------- +function out = node_root(xs, C) +xs_n = xs(:,C); +R = xs_n*xs_n'; +m = sum(xs_n,2); +N = length(C); + +node = struct('R',R,'m',m,'N',N,'C',C,'q',0,'e',0,'lambda',0, 'tse', 0); + +out = node_stat(node, xs); + +% --------------------------------------------------------- +function out = node_stat(in, xs) +xs_n = xs(:,in.C); + +Cov_n = in.R - 1/in.N*in.m*in.m'; +[ev,es] = eig(Cov_n); +[max_v,max_i] = max(diag(es)); +e_n = ev(:,max_i); +q_n = in.m/in.N; + +tse_n = norm(xs_n-repmat(q_n,1,in.N))^2; +lambda_n = sum(((xs_n - repmat(q_n,1,in.N))'*e_n).^2); +in.e = e_n; +in.q = q_n; +in.tse = tse_n; +in.lambda = lambda_n; +out = in; + + diff --git a/common/clusters_ransac.m b/common/clusters_ransac.m new file mode 100755 index 0000000..106e593 --- /dev/null +++ b/common/clusters_ransac.m @@ -0,0 +1,37 @@ +% clusters.m +% function cb = clusters(input,ref,maxIt) +% clusters(randn(2,400)/sqrt(12),randn(2,1)/sqrt(12)); +% 17.01.2003, Jens Ahrensfeld + +function [cb, idx, output] =clusters_ransac(input,max_distance,move_thresh) +cl_count = 0; +[nDim, nData] = size(input); +points_remain = input; + +while length(points_remain) > 0, + + [nDim,nRem] = size(points_remain); + new_idx = round(rand(1,1)*(nRem-1))+1; + ref = points_remain(:,new_idx); + + ref_last = zeros(nDim,1); + d_max = 2*move_thresh; + while d_max > move_thresh + d = sqrt(sum((points_remain-repmat(ref,1,nRem)).^2)); % Calc distances between references and input + points_assigned_idx = find(d <= max_distance); + num_assigned = length(points_assigned_idx); + if (num_assigned > 0); + ref = mean(points_remain(:,points_assigned_idx),2); + d_max = max(sqrt(sum((ref_last-ref).^2))); + ref_last = ref; + else + break; + end; + end; + cl_count = cl_count +1; + output{cl_count} = points_remain(:,points_assigned_idx); + idx{cl_count} = points_assigned_idx; + points_remain(:,points_assigned_idx) = []; + cb(:,cl_count) = ref; + +end; diff --git a/common/coef.m b/common/coef.m new file mode 100755 index 0000000..8fd5ac1 --- /dev/null +++ b/common/coef.m @@ -0,0 +1,71 @@ +% COEF generiert IIR-Tiefpass Koeffizienten für geradzahlige Ordnung +% [b,a] = coef(fa,fg,Q,N) +% +% fa : Abtastfrequenz [Hz] +% fg : Grenzfrequenz [Hz] fg <= fa/2 +% Q : Polgüte (Q=1.0 für Butterworth) +% N : Filterordnung (durch zwei teilbar) +% Type : 'Lowpass', 'Highpass', 'Bandpass', 'Bandstop' +% +function [b,a] = coef(fa,fg,Q,N,Type) + +if (rem(N,2) > 0) + error('Filterordnung muss z.Zt. noch durch zwei teilbar sein!'); +end; + +if (strcmp(Type,'Bandstop')) + TF = 4; +end; + +TF +sn = sin(2*pi*fg/fa); +cs = cos(2*pi*fg/fa); +ak = ones(N/2,3); +bk = ones(N/2,3); +a = 1; +b = 1; + +if (strcmp(Type,'Lowpass')) + for n = 1:N/2, + Qp = 1/(2*sin(pi*(n-0.5)/N)) + alpha = sn/(2*Q*Qp); + a0 = 1+alpha; + ak(n,2) = -2*cs/a0; + ak(n,3) = (1-alpha)/a0; + bk(n,1) = 0.5*(1-cs)/a0; + bk(n,2) = (1-cs)/a0; + bk(n,3) = 0.5*(1-cs)/a0; + a = conv(a,ak(n,(1:3))); + b = conv(b,bk(n,(1:3))); + end; +end; + +if (strcmp(Type,'Highpass')) + for n = 1:N/2, + Qp = 1/(2*sin(pi*(n-0.5)/N)) + alpha = sn/(2*Q*Qp); + a0 = 1+alpha; + ak(n,2) = -2*cs/a0; + ak(n,3) = (1-alpha)/a0; + bk(n,1) = 0.5*(1+cs)/a0; + bk(n,2) = -(1+cs)/a0; + bk(n,3) = 0.5*(1+cs)/a0; + a = conv(a,ak(n,(1:3))); + b = conv(b,bk(n,(1:3))); + end; +end; + +if (strcmp(Type,'Bandpass')) + for n = 1:N/2, + Qp = 1/(2*sin(pi*(n-0.5)/N)) + alpha = sn/(2*Q*Qp); + a0 = 1+alpha; + ak(n,2) = -2*cs/a0; + ak(n,3) = (1-alpha)/a0; + bk(n,1) = alpha/a0; + bk(n,2) = 0; + bk(n,3) = -alpha/a0; + a = conv(a,ak(n,(1:3))); + b = conv(b,bk(n,(1:3))); + end; +end; diff --git a/common/color_conv.m b/common/color_conv.m new file mode 100755 index 0000000..6065ed8 --- /dev/null +++ b/common/color_conv.m @@ -0,0 +1,8 @@ +function out = color_conv(in, type_in, type_out) +[ny,nx,nc] = size(in) + +% Get LAB-Data +color_write('color_in.dat',reshape(in,nc,ny*nx)); +commandStr = sprintf('color_conv.exe %s %s %s %s','color_in.dat', 'color_out.dat',type_in, type_out); +dos(commandStr); +out = reshape(color_read('color_out.dat'), ny,nx,nc); diff --git a/common/color_conv_v.m b/common/color_conv_v.m new file mode 100755 index 0000000..d2c88d7 --- /dev/null +++ b/common/color_conv_v.m @@ -0,0 +1,8 @@ +function out = color_conv(in, type_in, type_out) +[ny,nx] = size(in) + +% Get LAB-Data +color_write('color_in.dat',in); +commandStr = sprintf('color_conv.exe %s %s %s %s','color_in.dat', 'color_out.dat',type_in, type_out); +dos(commandStr); +out = color_read('color_out.dat'); diff --git a/common/color_read.m b/common/color_read.m new file mode 100755 index 0000000..9182f95 --- /dev/null +++ b/common/color_read.m @@ -0,0 +1,20 @@ +function data = color_read(name) +item_size = 8; + +fid = fopen(name,'r'); +if(fid < 0) + data = fid; + return; +end; + +version = fread(fid,1,'uint'); +prec_type = fscanf(fid,'%s',1); +fread(fid,1,'uint8'); +n_rows = fread(fid,1,'uint'); +n_cols = fread(fid,1,'uint'); +data_size = fread(fid,1,'uint'); +data = reshape(fread(fid,n_rows*n_cols,prec_type),n_rows,n_cols); +[n_rows,n_cols] = size(data) + +fclose(fid); + diff --git a/common/color_write.m b/common/color_write.m new file mode 100755 index 0000000..e7f354b --- /dev/null +++ b/common/color_write.m @@ -0,0 +1,16 @@ +function color_write(name, data) +[n_rows,n_cols] = size(data) + +item_size = 8; +prec_type = 'float64'; + +fid = fopen(name,'w'); +fwrite(fid,1,'uint'); +fprintf(fid,'%s',prec_type); +fwrite(fid,13,'uint8'); +fwrite(fid,n_rows,'uint'); +fwrite(fid,n_cols,'uint'); +fwrite(fid,n_rows*n_cols*item_size,'uint'); +fwrite(fid,data(:),prec_type); +fclose(fid); + diff --git a/common/correlt.m b/common/correlt.m new file mode 100755 index 0000000..33ee2ac --- /dev/null +++ b/common/correlt.m @@ -0,0 +1,15 @@ +% correl(x,y,N) + +function rxy = correlt(x,y,S) + +L = lge(x); +N = L-S; + +if (N == 0) + error('N = 0 !'); +end; + +for s=1:S, +rxy(s)=sum(x(S-s+1:L-s+1).*y(S:L))^2/(sum(x(S-s+1:L-s+1).^2)*sum(y(S:L).^2)); +%rxy(s)=abs(sum(x(S-s+1:L-s+1).*y(S:L)))/(sum(abs(x(S-s+1:L-s+1).*y(S:L)))+0.0001); +end; diff --git a/common/corrja.m b/common/corrja.m new file mode 100755 index 0000000..9aa9a81 --- /dev/null +++ b/common/corrja.m @@ -0,0 +1,15 @@ +% correl(x,y,N) + +function rxy = corrph(x,y,S) + +L = lge(x); +N = L-S; + +if (N == 0) + error('N = 0 !'); +end; + +for s=1:S, +%rxy(s)=sum(x(S-s+1:L-s+1).*y(S:L))^2/(sum(x(S-s+1:L-s+1).^2)*sum(y(S:L).^2)); +rxy(s)=abs(sum(x(S-s+1:L-s+1).*y(S:L)))/(sum(abs(x(1:L).*y(1:L)))+0.0001); +end; diff --git a/common/corrlp.m b/common/corrlp.m new file mode 100755 index 0000000..8af1e41 --- /dev/null +++ b/common/corrlp.m @@ -0,0 +1,15 @@ +% correl(x,y,N) + +function rxy = corrlp(x,y,S) + +L = lge(x); +N = L-S; + +if (N == 0) + error('N = 0 !'); +end; + +for s=1:S, +rxy(s)=sum(x(S-s+1:L-s+1).*y(S:L))^2/(sum(x(S-s+1:L-s+1).^2)*sum(y(S:L).^2)); +%rxy(s)=abs(sum(x(S-s+1:L-s+1).*y(S:L)))/(sum(abs(x(S-s+1:L-s+1).*y(S:L)))+0.0001); +end; diff --git a/common/corrph.m b/common/corrph.m new file mode 100755 index 0000000..737850b --- /dev/null +++ b/common/corrph.m @@ -0,0 +1,17 @@ +% correl(x,y,N) + +function rxy = corrph(x,y,S) + +L = lge(x); +N = L-S; + +if (N == 0) + error('N = 0 !'); +end; + +for s=0:S-1, +%rxy(s)=sum(x(S-s+1:L-s+1).*y(S:L))^2/(sum(x(S-s+1:L-s+1).^2)*sum(y(S:L).^2)); +%rxy(s)=abs(sum(x(S-s+1:L-s+1).*y(S:L)))/(sum(abs(x(S-s+1:L-s+1).*y(S:L)))+0.0001); +rxy(s+1)=sum(x(1:L-S).*y(1+s:L-S+s))/(sum(abs(x(1:L-S).*y(1+s:L-S+s)))+0.001); +%rxy(s+1)=sum((x(1:L-S).^2).*(y(1+s:L-S+s).^2))/(sum(x(1:L-S).^2 + y(1+s:L-S+s).^2) +0.001); +end; diff --git a/common/dat2wav.m b/common/dat2wav.m new file mode 100755 index 0000000..83391b2 --- /dev/null +++ b/common/dat2wav.m @@ -0,0 +1,17 @@ +% function y = dat2wav(name, fa, nBits, ampthresh) +% + +function y = dat2wav(name, fa, nBits, ampthresh) + +datfile = sprintf('%s.dat', name); +wavfile = sprintf('%s.wav', name); + +fid = fopen(datfile,'r'); +data_float = fread(fid, 'float32'); +fclose(fid); + +a = max(abs(data_float)); +ks = ampthresh/a; + +data = (2^(nBits-1))*ks*data_float; +wavw16(data,fa,wavfile); diff --git a/common/dosa.m b/common/dosa.m new file mode 100755 index 0000000..f8595a3 --- /dev/null +++ b/common/dosa.m @@ -0,0 +1,10 @@ +% correl(x,y,N) + +function xm = dosa(x,fa,M) + +xf = iirtp(fa,0.45*fa/M,1,8,x); + +xm = xf(1:M:lge(x)); + + + diff --git a/common/dspfft.m b/common/dspfft.m new file mode 100755 index 0000000..9140ad1 --- /dev/null +++ b/common/dspfft.m @@ -0,0 +1,3 @@ +% X = DSPFFT(x,C) +function X = dspfft(x,C) +X = fft(x,C)/C; diff --git a/common/dspifft.m b/common/dspifft.m new file mode 100755 index 0000000..367f3ac --- /dev/null +++ b/common/dspifft.m @@ -0,0 +1,3 @@ +% x = DSPIFFT(X,C) +function x = dspifft(X,C) +x = ifft(X,C)*C; diff --git a/common/dtdet01.m b/common/dtdet01.m new file mode 100755 index 0000000..da6aa8b --- /dev/null +++ b/common/dtdet01.m @@ -0,0 +1,85 @@ +% dtdet3 +function [a,rxds,xs,xl,es,clm,dtf] = dtdet(x,d,e,dti) +nstates = 9; +xsi = 1; +xli = 2; +dtsi = 3; +esi = 4; +nli = 5; +eli = 6; +rxdi = 7; +psxi = 8; +psdi = 9; +dtalki = 10; + +% ------------------------------------------------------------------------ +% Double-Talk-Detector params +ar = 0.05; +af = 0.005; +arrxd = 0.0004; +afrxd = 1.0; +al = 0.0001; +pa = 0.7; +pl = 2.0; +po = 0.93; +htime = 80; + +% ------------------------------------------------------------------------ +if (lge(dti) == 1); + dti = zeros(nstates,1); + dti(xli)=0.01; + dti(nli)=0.01; + dti(eli)=0.1; + dti(dtalki) = 0; +end; +% ------------------------------------------------------------------------ + +fa = 8000; +fc = 800; +L = lge(x); +M = 4; +% ------------------------------------------------------------------------ + +%d = iirtp(fa,fc,1,8,d); +%x = iirtp(fa,fc,1,8,x); + +d = dosa(d,fa,M); +x = dosa(x,fa,M); +e = dosa(e,fa,M); + + +%[xs,dti(xsi)] = phxs(x,ar,af,dti(xsi)); +%[xl,dti(xli)] = phxl(x,xs,al,pl,dti(xli)); + + + +rxd = max(abs(corrph(x,d,40/M))); + +%[es,dti(esi)] = phxs(e,0.008,0.005,dti(esi)); +%[el,dti(eli)] = phxl(e,es,0.0001,pl,dti(eli)); +%[rxds,dti(rxdi)] = phxs(rxd*ones(L,1),arrxd,afrxd,dti(rxdi)); + +%[nl,dti(nli)] = phnl(e,es,el,xs,xl,al,pl,dti(nli)); +%a = alpha(xs,xl,nl,pl,es,(rxd > po),pa); + + +if (rxd > po) + if(dti(dtalki)==1) + a = [zeros(htime,1);ones(L-htime,1)]; + dti(dtalki)=0; + else + a = ones(L,1); + end; +else + a = zeros(L,1); + dti(dtalki)=1; +end; +rxds=ones(L,1)*rxd; + +dtf = dti; +xs=1; +es=1; +xl=1; +clm=1; + +% ------------------------------------------------------------------------ diff --git a/common/dtdet02.m b/common/dtdet02.m new file mode 100755 index 0000000..e2be6ab --- /dev/null +++ b/common/dtdet02.m @@ -0,0 +1,82 @@ +% dtdet3 +function [a,rxds,xs,xl,ds,dl,dtf] = dtdet(x,d,e,dti) +nstates = 11; +xsi = 1; +xli = 2; +dtsi = 3; +dsi = 4; +dli = 5; +eli = 6; +rxdi = 7; +psxi = 8; +psdi = 9; +dtalki = 10; +esi = 11; +% ------------------------------------------------------------------------ +% Double-Talk-Detector params +ar = 0.05; +af = 0.005; +arrxd = 0.0004; +afrxd = 1.0; +al = 0.0001; +pa = 0.7; +pl = 2.0; +po = 0.93; +arf = 0.0006; +htime = 80; + +% ------------------------------------------------------------------------ +if (lge(dti) == 1); + dti = zeros(nstates,1); + dti(xli)=0.1; + dti(dli)=0.1; + dti(xsi)=0.0; + dti(dsi)=0.0; + dti(eli)=0.1; + dti(dtalki) = 0; +end; +% ------------------------------------------------------------------------ + +fa = 8000; +fc = 800; +L = lge(x); +M = 1; +% ------------------------------------------------------------------------ + +%d = iirtp(fa,fc,1,8,d); +%x = iirtp(fa,fc,1,8,x); + +d = dosa(d,fa,M); +x = dosa(x,fa,M); +e = dosa(e,fa,M); + + +[xs,dti(xsi)] = sm2(x,dti(xli),arf); +[xl,dti(xli)] = sm2(x,dti(xli),arf/M); +%[xl,dti(xli)] = phxl(x,xs,al,pl,dti(xli)); +%[ds,dti(dsi)] = phxs(d,ar,af,dti(dsi)); +[dl,dti(dli)] = sm2(d,dti(dli),arf/M); +[el,dti(eli)] = sm2(e,dti(eli),0.5*arf/M); +%[es,dti(esi)] = sm2(e,dti(esi),10*arf/M); +%[dl,dti(dli)] = phxl(d,ds,al,pl,dti(dli)); +xs = zeros(L,1); +ds = zeros(L,1); + +%rxd = max(abs(correlt(x,d,40/M))); + +%[es,dti(esi)] = phxs(e,0.008,0.005,dti(esi)); +%[el,dti(eli)] = phxl(e,es,0.0001,pl,dti(eli)); +%[rxds,dti(rxdi)] = phxs(rxd*ones(L,1),arrxd,afrxd,dti(rxdi)); + +%[nl,dti(nli)] = phnl(e,es,el,xs,xl,al,pl,dti(nli)); +%a = alpha(xs,xl,nl,pl,es,(rxd > po),pa); +a = zeros(L,1); + +rxd = max(corrph(xl,dl,30)); +a = (((xl(1:1.0/M:L/M)*rxd) > (pl*dl(1:1.0/M:L/M)))); +rxds=0*ones(L,1); +dtf = dti; +xl = zeros(L,1); +dl = zeros(L,1); + +% ------------------------------------------------------------------------ diff --git a/common/dtdet03.m b/common/dtdet03.m new file mode 100755 index 0000000..b1280ac --- /dev/null +++ b/common/dtdet03.m @@ -0,0 +1,85 @@ +% dtdet3 +function [a,rxds,xs,xl,es,clm,dtf] = dtdet(x,d,e,dti) +nstates = 9; +xsi = 1; +xli = 2; +dtsi = 3; +esi = 4; +nli = 5; +eli = 6; +rxdi = 7; +psxi = 8; +psdi = 9; +dtalki = 10; + +% ------------------------------------------------------------------------ +% Double-Talk-Detector params +ar = 0.05; +af = 0.005; +arrxd = 0.0004; +afrxd = 1.0; +al = 0.0001; +pa = 0.7; +pl = 2.0; +po = 0.93; +htime = 80; + +% ------------------------------------------------------------------------ +if (lge(dti) == 1); + dti = zeros(nstates,1); + dti(xli)=0.01; + dti(nli)=0.01; + dti(eli)=0.1; + dti(dtalki) = 0; +end; +% ------------------------------------------------------------------------ + +fa = 8000; +fc = 800; +L = lge(x); +M = 4; +% ------------------------------------------------------------------------ + +%d = iirtp(fa,fc,1,8,d); +%x = iirtp(fa,fc,1,8,x); + +d = dosa(d,fa,M); +x = dosa(x,fa,M); +e = dosa(e,fa,M); + + +%[xs,dti(xsi)] = phxs(x,ar,af,dti(xsi)); +%[xl,dti(xli)] = phxl(x,xs,al,pl,dti(xli)); + + + +rxd = min(max(corrja(x,d,30)),1); + +%[es,dti(esi)] = phxs(e,0.008,0.005,dti(esi)); +%[el,dti(eli)] = phxl(e,es,0.0001,pl,dti(eli)); +%[rxds,dti(rxdi)] = phxs(rxd*ones(L,1),arrxd,afrxd,dti(rxdi)); + +%[nl,dti(nli)] = phnl(e,es,el,xs,xl,al,pl,dti(nli)); +%a = alpha(xs,xl,nl,pl,es,(rxd > po),pa); + + +if (rxd > po) + if(dti(dtalki)==1) + a = [zeros(htime,1);ones(L-htime,1)]; + dti(dtalki)=0; + else + a = ones(L,1); + end; +else + a = zeros(L,1); + dti(dtalki)=1; +end; +rxds=ones(L,1)*rxd; + +dtf = dti; +xs=1; +es=1; +xl=1; +clm=1; + +% ------------------------------------------------------------------------ diff --git a/common/eudist.m b/common/eudist.m new file mode 100755 index 0000000..4fe2b95 --- /dev/null +++ b/common/eudist.m @@ -0,0 +1,31 @@ +% clusters.m +% function cb = clusters(input,ref,maxIt) +% clusters(randn(2,400)/sqrt(12),randn(2,1)/sqrt(12)); +% 17.01.2003, Jens Ahrensfeld + +function [indices, output] = eudist(input,ref) + +[nDim, nData] = size(input); +[nDimr, nClasses] = size(ref); + +if (nDim ~= nDimr), + error('Dimensions must match!'); +end; + +cnt = zeros(nClasses,1); +for j=1:nClasses, + indices{j} = []; +end; +for i=1:nData, + for j=1:nClasses, + d(j) = sqrt(sum((input(:,i)-ref(:,j)).^2)); % Calc distances between references and input + end; + [value nearestClass] = min(d); + cnt(nearestClass) = cnt(nearestClass) + 1; + indices{nearestClass} = [ indices{nearestClass} i ]; +end; + +% Separate data +for j=1:nClasses, + output{j} = input(:,indices{j}); +end; \ No newline at end of file diff --git a/common/flms.m b/common/flms.m new file mode 100755 index 0000000..d0ad940 --- /dev/null +++ b/common/flms.m @@ -0,0 +1,133 @@ +% Frequency LMS Adaptive Filter -FLMS- +% ------------------------------------ +% +% Aufruf: +% Function [ys,ee,wf,fpo]=flms(x,d,alpha,N,C,L,wi) +% +% Parameter: +% x : Eingangsvektor +% d : Referenzvektor +% alpha: konstante Schrittweite 0 < alpha <= 1.0 +% N : Anzahl Filterkoeffizienten +% C : FFT-Länge, idealerweise eine Zweier-Potenz (0=Auto) +% L : Länge eines Verarbeitungsblocks (idealerweise C=L+N=2N) +% wi : Startwerte der Filterkoeffizienten im Zeitbereich +% +% Rückgabewerte: +% ee : Adaptitionsfehler +% ys : Filterausgang +% wf : Filterkoeffizienten am Ende der Adaption im Zeitbereich +% fpo : Anzahl benötigten der Floating-Point Operationen + +% ------------------------------------------------------------------------- +% Datum : 20.7.2001 +% Autor : Jens Ahrensfeld +% Thema : Diplomarbeit +% Datei : flms.m +% Benötigte Dateien: dspfft.m, dspifft.m +% +% ------------------------------------------------------------------------- +function [ys,ee,wf,fpo]=flms(x,d,alpha,N,C,L,wi) + +% ------------------------------------------------------------------------- +% Überprüfen der Parameter +% ------------------------------------------------------------------------- +if (N ~= fix(N)), + error('Fehler: N muss eine ganze Zahl sein!'); +end; + +if (C ~= fix(C)), + error('Fehler: C muss eine ganze Zahl sein!'); +end; + +if (L ~= fix(L)), + error('Fehler: L muss eine ganze Zahl sein!'); +end; + +% ------------------------------------------------------------------------- +% Automatik-Auswahl +% ------------------------------------------------------------------------- +% Falls C=0 => Automatische Wahl von C=2^r +% Anahme: C = L + N +if (C==0) + C = 2^ceil(log2(L+N-1)) +end; + +% ------------------------------------------------------------------------- +% FLMS Parameter +% ------------------------------------------------------------------------- +gamma = 1.0/C; % Sicherheitskonstante vermeidet Division durch Null +lambda = 0.6; % Vergessensfaktor + +% ------------------------------------------------------------------------- +% Initialisierungen +% ------------------------------------------------------------------------- +Lx = length(x); % Länge des Eingangsvektors +K = Lx/L; % Anzahl Verarbeitungsblöcke + +ee=zeros(Lx,1); % Fehlervektor +ys=zeros(Lx,1); % Filterausgang, Schaetzung von d +wf = zeros(C,1); % Filter im Zeitbereich +PX = zeros(C,1); % Schaetzung des Leistungsdichtespektrums + +WS = zeros(C,1); % Filterkoeff. im Frequenzbereich +X = zeros(C,1); % Eingangsvektor im Frequenzbereich +mu = zeros(C,1); % Variable Schrittweite +yy = zeros(C,1); % Filterausgang im Zeitbereich +YS = zeros(C,1); % Filterausgang im Frequenzbereich +xzp = zeros(Lx+C-L,1); % Eingangsvektor mit C-L führenden Nullen + +xzp(C-L+1:Lx+C-L) = x; % Einfügen von Np führenden Nullen an den Eingangsvektor + +% Anfangswerte der Filterkoeffienten in den Frequenzbereich transformieren +WS(1:C) = dspfft(wi(1:N),C); + +% ------------------------------------------------------------------------- +% Algorithmus Start +% ------------------------------------------------------------------------- +flops(0); +for k=1:K, + + kL = (k-1)*L; + + % Transformation des aktuellen Eingangsvektors in den Frequenzbereich + X(1:C) = dspfft(xzp(kL+1:kL+C),C); + + % Schnelle FFT-Faltung + YS = WS(1:C).*X(1:C) *C; + + % Transformation des Filterausgangs in den Zeitbereich + yy = real(dspifft(YS,C)); + + % Abspeichern der letzten L Werte des Filterausgangs + ys(kL+1:kL+L) = yy(C-L+1:C); + + % Berechnung des Fehlers e(k) + ee(kL+1:kL+L) = d(kL+1:kL+L) - ys(kL+1:kL+L); + + % Transformation des Fehlers in den Frequenzbereich + E = dspfft([zeros(C-L,1); ee(kL+1:kL+L)],C); + + % Schätzung der mittleren Eingangsleistung Px(k) aus x(k) + PX = abs((1-lambda)*conj(X(1:C)).*X(1:C) *C + lambda*PX); + + % Berechnung der variablen Schrittweite mu(k) + % mu(k) wird niemals grösser als eins + mu(1:C) = (alpha*gamma) ./(PX+gamma); + + % Filterkoeffizenten-Update + WS(1:C) = WS(1:C) + mu(1:C) .* conj(X(1:C)) .* E *C; + + +% Projektion der Filterkoeffizienten + wf(1:C) = real(dspifft(WS(1:C),C)); + WS(1:C) = dspfft(wf(1:N),C); + +end; +wf = [wf(1:C-L+1) ; zeros(L-1,1)]; % Filterkoeffizienten am Ende der Adaption + +fpo=flops; + +% ------------------------------------------------------------------------- +% Ende FLMS.M +% ------------------------------------------------------------------------- diff --git a/common/flms1.m b/common/flms1.m new file mode 100755 index 0000000..3497bc7 --- /dev/null +++ b/common/flms1.m @@ -0,0 +1,60 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% Fast LMS Algorithm % +% % +% Written By: Sundar Sankaran and A. A. (Louis) Beex % +% DSP Research Laboratory % +% Dept. of Electrical and Comp. Engg % +% Virginia Tech % +% Blacksburg VA 24061-0111 % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +randn('seed', 0) ; +rand('seed', 0) ; + +NoOfData = 8000 ; % Set no of data points used for training +M = 32 ; % Set the adaptive filter order + +Mu = 0.01 ; % Set the step-size constant +Gamma = 0.9 ; % Forgetting factor +Delta = 0.01 ; % R_est initialized to Delta*I + +u = randn(NoOfData, 1) ;% Input assumed to be white +h = rand(M, 1) ; % System picked randomly +d = filter(h, 1, u) ; % Generate output (desired signal) + +% Initialize fast-lms + +W = zeros(2*M,1) ; +p = Delta*ones(2*M,1) ; + +y = zeros(M, 1) ; +e = zeros(M, 1) ; + +for k = 2 : floor(length(u)/M) - 1 ; + + U = fft(u((k-1)*M:(k+1)*M-1)) ; + Y = ifft(U.*W) ; + y(k*M:(k+1)*M-1) = real(Y(M+1:2*M)) ; + + e(k*M:(k+1)*M-1) = d(k*M:(k+1)*M-1) - y(k*M:(k+1)*M-1) ; + E = fft([zeros(M,1); e(k*M:(k+1)*M-1)]) ; + + p = Gamma*p+(1-Gamma)*abs(U).^2 ; + p_inv = 1./p ; + + PHI = ifft(p .* conj(U) .* E) ; + phi = real(PHI(1 : M )) ; + W = W + Mu / (2*M) * fft([phi; zeros(M,1)]) ; + +end ; + +% Plot results + +figure ; +plot(20*log10(abs(e))) ; +title('Learning Curve') ; +xlabel('Iteration Number') ; +ylabel('Output Estimation Error in dB') ; + diff --git a/common/flms2.m b/common/flms2.m new file mode 100755 index 0000000..f630b50 --- /dev/null +++ b/common/flms2.m @@ -0,0 +1,99 @@ +%=================================================================== +%------------------------------------------------------------------- +% +% Adaptive Filter +% +% Simulationen +% +% [ee,w,dd,calc,Px]=flms2(mu,N,C,L,x,d,w_start) +% +% FLMS-Algorithmus +% +% ee = Fehlersignal +% w = Filterkoeffizienten am Ende der Adaption +% dd = Filterausgang, Schaetzung von d +% calc = Rechenzeit +% Px = Schaetzung des Leistungsdichtespektrums vom x +% mu = Schrittweite +% N = Anz. Filterkoeffizienten = Filterordnung+1 +% X = Filtereingang x[.] +% d = erwuenschtes Signal d[.] +% w_start = Startwerte Filterkoeffizienten +% +%------------------------------------------------------------------- + +% +% author: Markus Hofbauer +% ISI, ETH Zuerich (Switzerland) +% +% created: 7/2000 +% +% +%------------------------------------------------------------------- +% +% File : flms2.m +% +% Startfile: sim5.m +% +%------------------------------------------------------------------- +%=================================================================== + +function [ee,w,dd,calc,PX]=flms2(mu,N,C,L,x,d,w_start) + +% Initialisierungen + +ee=zeros(lge(x),1); % Fehlervektor +dd=zeros(lge(x),1); % Filterausgang, Schaetzung von d +WS = zeros(C,1); % Gewichtsvektor im Frequenzbereich +w = zeros(C,1); % Gewichtsvektor im Zeitbereich +PX = C*ones(C,1); % Schaetzung des Leistungsdichtespektrums +X = zeros(C,1); % Eingangsvektor im Frequenzbereich + +%% FLMS Parameter +numblk = floor(lge(x)/L); % Anzahl Verarbeitungsbloecke +begin=ceil(C/L)-1; % erster Blockindex +gamma = 0.6; % Vergessensfaktor + +tic; % Stopuhr laeuft + +%------------------------------------------------------------------- +% FLMS Update Loop +%------------------------------------------------------------------- + +for j = begin: numblk-1 + + % DFT von x, Laenge C + X = fft(x((j+1)*L-C+1:(j+1)*L)); % aktuellster Wert + % des Blockes: (j+1)*L + % Filterung im Frequenzbereich + YS = (X .* WS); + ys = real(ifft(YS)); % IDFT, Filterausgang im Zeitbereich + + % Fehlersignal im Zeitbereich + ee(j*L+1:(j+1)*L) = d(j*L+1:(j+1)*L) - ys(C-L+1:C); + % Overlap-Save-Verfahren: nur die letzten L Werte von ys + % werden verwendet um eine lineare Faltung zu erhalten + + % Fehlersignal im Frequenzbereich + E = fft([zeros(C-L,1); ee(j*L+1:(j+1)*L)]); + + % Update von PX + PX=abs((1-gamma)*conj(X).*X+ gamma * PX); + + % Update von WS + WS = WS + mu * conj(X)./(PX+0.001) .* E ; % Adaption von WS + + % Projektion (letzten L-1 Werte von w zu Null setzten) + w = real(ifft(WS)); + WS = fft(w(1:C-L+1),C); + + + dd(j*L+1:(j+1)*L) = ys(C-L+1:C); % Abspeichern der Schaetzung + +end; + +calc=toc; % Stopuhr angehalten +w = [w(1:C-L+1) ; zeros(L-1,1)]; % Filterkoeffizienten am Ende der Adaption + + + diff --git a/common/flms3.m b/common/flms3.m new file mode 100755 index 0000000..cf5fceb --- /dev/null +++ b/common/flms3.m @@ -0,0 +1,108 @@ +%=================================================================== +%------------------------------------------------------------------- +% +% Adaptive Filter +% +% Simulationen +% +% [ee,w,dd,calc,Px]=flms2(mu,N,C,L,x,d,w_start) +% +% FLMS-Algorithmus +% +% ee = Fehlersignal +% w = Filterkoeffizienten am Ende der Adaption +% dd = Filterausgang, Schaetzung von d +% calc = Rechenzeit +% Px = Schaetzung des Leistungsdichtespektrums vom x +% mu = Schrittweite +% N = Anz. Filterkoeffizienten = Filterordnung+1 +% X = Filtereingang x[.] +% d = erwuenschtes Signal d[.] +% w_start = Startwerte Filterkoeffizienten +% +%------------------------------------------------------------------- + +% +% author: Markus Hofbauer +% ISI, ETH Zuerich (Switzerland) +% +% created: 7/2000 +% +% +%------------------------------------------------------------------- +% +% File : flms2.m +% +% Startfile: sim5.m +% +%------------------------------------------------------------------- +%=================================================================== + +function [ee,w,dd,calc,PX,u]=flms2(mu,N,C,L,x,d,w_start) + +% Initialisierungen + +ee=zeros(lge(x),1); % Fehlervektor +dd=zeros(lge(x),1); % Filterausgang, Schaetzung von d +WS = zeros(C,1); % Gewichtsvektor im Frequenzbereich +w = zeros(C,1); % Gewichtsvektor im Zeitbereich +PX = zeros(C,1); % Schaetzung des Leistungsdichtespektrums +X = zeros(C,1); % Eingangsvektor im Frequenzbereich +u = zeros(C,1); + +%% FLMS Parameter +numblk = floor(lge(x)/L); % Anzahl Verarbeitungsbloecke +begin=ceil(C/L)-1; % erster Blockindex +gamma = 0.97; % Vergessensfaktor +alpha = 1.0; % 0 < alpha < 1 + +tic; % Stopuhr laeuft + +dspMul = C +dspShift = log2(dspMul) + +%------------------------------------------------------------------- +% FLMS Update Loop +%------------------------------------------------------------------- + +for j = begin: numblk-1 + + % DFT von x, Laenge C + X = dspfft(x((j+1)*L-C+1:(j+1)*L),C); % aktuellster Wert + % des Blockes: (j+1)*L + % Filterung im Frequenzbereich + YS = ((X*dspMul) .* (WS*dspMul)); + + ys = real(dspifft(YS,C)); % IDFT, Filterausgang im Zeitbereich + + % Fehlersignal im Zeitbereich + ee(j*L+1:(j+1)*L) = d(j*L+1:(j+1)*L) - ys(C-L+1:C); + % Overlap-Save-Verfahren: nur die letzten L Werte von ys + % werden verwendet um eine lineare Faltung zu erhalten + + % Fehlersignal im Frequenzbereich + E = dspfft([zeros(C-L,1); ee(j*L+1:(j+1)*L)],C); + + % Update von PX + PX=abs((1-gamma)*conj(X).*X*dspMul + gamma * PX); + + % umax = alpha + u = (alpha*mu/dspMul) ./(PX+mu/dspMul); + + % Update von WS + WS = WS + u .* conj(X).* E ; % Adaption von WS + + % Projektion (letzten L-1 Werte von w zu Null setzten) + w = real(dspifft(WS,C)); + WS = dspfft(w(1:C-L+1),C); + + + dd(j*L+1:(j+1)*L) = ys(C-L+1:C); % Abspeichern der Schaetzung + +end; + +calc=toc; % Stopuhr angehalten +w = [w(1:C-L+1) ; zeros(L-1,1)]; % Filterkoeffizienten am Ende der Adaption + + + diff --git a/common/flms_d.m b/common/flms_d.m new file mode 100755 index 0000000..51e913d --- /dev/null +++ b/common/flms_d.m @@ -0,0 +1,142 @@ +% Frequency LMS Adaptive Filter -FLMS- (DEBUG) +% -------------------------------------------- +% +% Aufruf: +% Function [ys,ee,wf,fpo]=flms_d(x,d,alpha,N,C,L,wi) +% +% Parameter: +% x : Eingangsvektor +% d : Referenzvektor +% alpha: konstante Schrittweite 0 < alpha <= 1.0 +% N : Anzahl Filterkoeffizienten +% C : FFT-Länge, idealerweise eine Zweier-Potenz (0=Auto) +% L : Länge eines Verarbeitungsblocks (idealerweise C=L+N=2N) +% wi : Startwerte der Filterkoeffizienten im Zeitbereich +% +% Rückgabewerte: +% ee : Adaptitionsfehler +% ys : Filterausgang +% wf : Filterkoeffizienten am Ende der Adaption +% fpo : Anzahl benötigten der Floating-Point Operationen +% +% Zu Debug-Zwecken werden folgende Variablen als Binär-Dateien gespeichert: +% wf.dat : Filterkoeffizienten am Ende der Adaption (Länge N) + +% ------------------------------------------------------------------------- +% Datum : 20.7.2001 +% Autor : Jens Ahrensfeld +% Thema : Diplomarbeit +% Datei : flms_d.m +% Benötigte Dateien: dspfft.m, dspifft.m +% +% ------------------------------------------------------------------------- +function [ys,ee,wf,fpo]=flms_d(x,d,alpha,N,C,L,wi) + +% ------------------------------------------------------------------------- +% Überprüfen der Parameter +% ------------------------------------------------------------------------- +if (N ~= fix(N)), + error('Fehler: N muss eine ganze Zahl sein!'); +end; + +if (C ~= fix(C)), + error('Fehler: C muss eine ganze Zahl sein!'); +end; + +if (L ~= fix(L)), + error('Fehler: L muss eine ganze Zahl sein!'); +end; + +% ------------------------------------------------------------------------- +% Automatik-Auswahl +% ------------------------------------------------------------------------- +% Falls C=0 => Automatische Wahl von C=2^r +% Anahme: C = L + N +if (C==0) + C = 2^ceil(log2(L+N-1)) +end; + +% ------------------------------------------------------------------------- +% FLMS Parameter +% ------------------------------------------------------------------------- +gamma = 1.0/C; % Sicherheitskonstante vermeidet Division durch Null +lambda = 0.6; % Vergessensfaktor + +% ------------------------------------------------------------------------- +% Initialisierungen +% ------------------------------------------------------------------------- +Lx = length(x); % Länge des Eingangsvektors +K = Lx/L; % Anzahl Verarbeitungsblöcke + +ee=zeros(Lx,1); % Fehlervektor +ys=zeros(Lx,1); % Filterausgang, Schaetzung von d +wf = zeros(C,1); % Filter im Zeitbereich +PX = zeros(C,1); % Schaetzung des Leistungsdichtespektrums + +WS = zeros(C,1); % Filterkoeff. im Frequenzbereich +X = zeros(C,1); % Eingangsvektor im Frequenzbereich +mu = zeros(C,1); % Variable Schrittweite +yy = zeros(C,1); % Filterausgang im Zeitbereich +YS = zeros(C,1); % Filterausgang im Frequenzbereich +xzp = zeros(Lx+C-L,1); % Eingangsvektor mit C-L führenden Nullen + +xzp(C-L+1:Lx+C-L) = x; % Einfügen von Np führenden Nullen an den Eingangsvektor + +% Anfangswerte der Filterkoeffienten in den Frequenzbereich transformieren +%WS(1:C) = dspfft(wi(1:N),C); + +% ------------------------------------------------------------------------- +% Algorithmus Start +% ------------------------------------------------------------------------- +flops(0); +for k=1:K, + + kL = (k-1)*L; + + % Transformation des aktuellen Eingangsvektors in den Frequenzbereich + X(1:C) = dspfft(xzp(kL+1:kL+C),C); + + % Schnelle FFT-Faltung + YS = WS(1:C).*X(1:C) *C; + + % Transformation des Filterausgangs in den Zeitbereich + yy = real(dspifft(YS,C)); + + % Abspeichern der letzten L Werte des Filterausgangs + ys(kL+1:kL+L) = yy(C-L+1:C); + + % Berechnung des Fehlers e(k) + ee(kL+1:kL+L) = d(kL+1:kL+L) - ys(kL+1:kL+L); + + % Transformation des Fehlers in den Frequenzbereich + E = dspfft([zeros(C-L,1); ee(kL+1:kL+L)],C); + + % Schätzung der mittleren Eingangsleistung Px(k) aus x(k) + PX = abs((1-lambda)*conj(X(1:C)).*X(1:C) *C + lambda*PX); + + % Berechnung der variablen Schrittweite mu(k) + % mu(k) wird niemals grösser als eins + mu(1:C) = (alpha*gamma) ./(PX+gamma); + + % Filterkoeffizenten-Update + WS(1:C) = WS(1:C) + mu(1:C) .* conj(X(1:C)) .* E *C; + + +% Projektion der Filterkoeffizienten + wf(1:C) = real(dspifft(WS(1:C),C)); + WS(1:C) = dspfft(wf(1:N),C); + +end; + +wf = [wf(1:C-L+1) ; zeros(L-1,1)]; % Filterkoeffizienten am Ende der Adaption + +% Speichern der Koeffizienten in Datei (Genauigkeit Float32) +fid = fopen('wf.dat','wb'); +fwrite(fid,wf,'float32'); +fclose(fid); + +fpo=flops; + +% ------------------------------------------------------------------------- +% Ende FLMS_D.M +% ------------------------------------------------------------------------- diff --git a/common/fpoints.m b/common/fpoints.m new file mode 100755 index 0000000..d03bc31 --- /dev/null +++ b/common/fpoints.m @@ -0,0 +1,11 @@ +% fpoints.m + +function fp = fpoints(dim,len, pos) + +fp = randn(dim,len); + +for d=1:dim, + fp(d,:) = fp(d,:) + pos(d); +end; + + diff --git a/common/iirtp.m b/common/iirtp.m new file mode 100755 index 0000000..47ee1ec --- /dev/null +++ b/common/iirtp.m @@ -0,0 +1,13 @@ +% IIR-Tiefpass +% y = iirtp(fa,fg,Q,N,x) +% +% fa : Abtastfrequenz [Hz] +% fg : Grenzfrequenz [Hz] fg <= fa/2 +% Q : Polgüte (Resonanz) +% N : Filterordnung (durch zwei teilbar) +% x : Eingangssignal +function y = iirtp(fa,fg,Q,N,x) + +[b,a] = iirtpc(fa,fg,Q,N); + +y = filter(b,a,x); diff --git a/common/iirtp2.m b/common/iirtp2.m new file mode 100755 index 0000000..d0d10da --- /dev/null +++ b/common/iirtp2.m @@ -0,0 +1,13 @@ +% IIR-Tiefpass +% [y,zf] = iirtp2(fa,fg,Q,N,x,zi) +% +% fa : Abtastfrequenz [Hz] +% fg : Grenzfrequenz [Hz] fg <= fa/2 +% Q : Polgüte (Resonanz) +% N : Filterordnung (durch zwei teilbar) +% x : Eingangssignal +function [y,zf] = iirtp(fa,fg,Q,N,x,zi) + +[b,a] = iirtpc(fa,fg,Q,N); + +[y,zf] = filter(b,a,x,zi); diff --git a/common/iirtpc.m b/common/iirtpc.m new file mode 100755 index 0000000..4d17c11 --- /dev/null +++ b/common/iirtpc.m @@ -0,0 +1,33 @@ +% IIRTPC generiert IIR-Tiefpass Koeffizienten für geradzahlige Ordnung +% [b,a] = iirtpc(fa,fg,Q,N) +% +% fa : Abtastfrequenz [Hz] +% fg : Grenzfrequenz [Hz] fg <= fa/2 +% Q : Polgüte (Resonanz) +% N : Filterordnung (durch zwei teilbar) +% +function [b,a] = iirtpc(fa,fg,Q,N) + +if (rem(N,2) > 0) + error('Filterordnung muss z.Zt. noch durch zwei teilbar sein!'); +end; + +sn = sin(2*pi*fg/fa); +cs = cos(2*pi*fg/fa); +ak = ones(N/2,3); +bk = ones(N/2,3); +a = 1; +b = 1; + +for n = 1:N/2, + Qp = 1/(2*sin(pi*(n-0.5)/N)); + alpha = sn/(2*Q*Qp); + a0 = 1+alpha; + ak(n,2) = -2*cs/a0; + ak(n,3) = (1-alpha)/a0; + bk(n,1) = 0.5*(1-cs)/a0; + bk(n,2) = (1-cs)/a0; + bk(n,3) = 0.5*(1-cs)/a0; + a = conv(a,ak(n,(1:3))); + b = conv(b,bk(n,(1:3))); +end; diff --git a/common/l1.m b/common/l1.m new file mode 100755 index 0000000..c80dd1b --- /dev/null +++ b/common/l1.m @@ -0,0 +1,48 @@ +% ################################################################################## +% ## Loesung: Impulsfolge ## +% ################################################################################## + +% ##### Teilaufgabe a: zeitdiskrete endliche Signale ##### +% ### 1: +k = 1:20; +x1 = zeros(size(k)); +x1(5) = 0.9; + +figure; stem(k,x1); grid; xlabel('k'); ylabel('x1(k)'); + +% ### 2: +k = -15:15; +x2 = zeros(size(k)); +x2(16) = 0.8; + +figure; stem(k,x2); grid; xlabel('k'); ylabel('x2(k)'); + +% ### 3: +k = 300:350; +x3 = zeros(size(k)); +x3(34) = 1.5; + +figure; stem(k,x3); grid; xlabel('k'); ylabel('x3(k)'); + +% ### 4: +k = -10:0; +x4 = zeros(size(k)); +x4(4) = 4.5; + +figure; stem(k,x4); grid; xlabel('k'); ylabel('x4(k)'); + +% ##### Teilaufgabe b: Impulsfolge ##### +P = 5; +M = 10; +k = 0:(M*P)-1; +x = [1; 0; 0; 0; 0] * ones(1,M); +s = x(:); + +figure; stem(k',s); grid; xlabel('k'); ylabel('s(k)'); axis([0 49 0 1.2]); + +% ##### Teilaufgabe c: Matlab Code ##### +x = [0;1;1;0;0;0] * ones(1,7); +x = x(:); + +figure; stem(x); xlabel('k'); ylabel('x(k)'); axis([0 42 0 1]); +% ##### EOF ##### diff --git a/common/l10.m b/common/l10.m new file mode 100755 index 0000000..81b588a --- /dev/null +++ b/common/l10.m @@ -0,0 +1,32 @@ +% ################################################################################## +% ## Loesung: Gleichverteilter zeitdiskreter stochastischer Prozess ## +% ################################################################################## +N=10000; + +% ##### Teilaufgabe a: Rauschfolge ##### +n = rand(N,1); + +% ##### Teilaufgabe b: Histogramm ##### +n = rand(N,1); +[H,x] = hist(n,30); +H = H/N; % Damit H eine Schaetzung der Wahrscheinlichkeitsdichte + % ist, muss die Flaeche eins sein + +figure; bar(x,H);grid; xlabel('x'); ylabel('f_X(x)'); +title(['Gleichvert. Mittelw.=',num2str(mean(n)),' Var.=',num2str(std(n).^2)]); + +% ##### Teilaufgabe c: Mittelwert und Varianz ##### +n = rand(N,1); +m = mean(n); +s = std(n); + +% ##### Teilaufgabe d: Musterfunktion ##### +for k=1:5 + n = rand(N,1); + [H,x]=hist(n,30); + H = H/N; + + figure; bar(x,H); grid; xlabel('x'); ylabel('f_X(x)'); + title(['Gleichvert. Mittelw.=',num2str(mean(n)),' Var.=',num2str(std(n).^2)]); +end; +% ##### EOF ##### diff --git a/common/l11.m b/common/l11.m new file mode 100755 index 0000000..39e7d10 --- /dev/null +++ b/common/l11.m @@ -0,0 +1,29 @@ +% ################################################################################## +% ## Loesung: Korrelationsfolgen ## +% ## -------------------------------------------------------------------------- ## +% ## Benoetigte(s) m-File(s): ldtft.m ## +% ################################################################################## + +% ##### Teilaufgabe a: weisse Rauschfolge ##### +n = randn(800,1); +kappa = -799:799; +s_biased = xcorr(n,'biased'); +s_unbiased = xcorr(n,'unbiased'); +s_coeff = xcorr(n,'coeff'); + +figure; stem(kappa,s_biased); grid; xlabel('kappa'); ylabel('Amplitude'); +title('Autokorrelationsfolge: Nicht erwartungstreu'); +figure; stem(kappa,s_unbiased); grid; xlabel('kappa'); ylabel('Amplitude'); +title('Autokorrelationsfolge: erwartungstreu'); +figure; stem(kappa,s_coeff); grid; xlabel('kappa'); ylabel('Amplitude'); +title('Autokorrelationsfolge: auf time-lag Null normiert'); + +% ##### Teilaufgabe b: Leistungsdichtespektrum ##### +n = randn(200,1); +s = xcorr(n,'biased'); +[S W] = ldtft(s,512); + +figure; semilogy(W/pi, abs(S) ); % normierte Frequenz +grid; xlabel('Omega/pi'); ylabel('|S_xx(exp(j*Omega))|'); +title('Leistungsdichtespektrum'); +% ##### EOF ##### diff --git a/common/l12.m b/common/l12.m new file mode 100755 index 0000000..6d04b27 --- /dev/null +++ b/common/l12.m @@ -0,0 +1,24 @@ +% ################################################################################## +% ## Loesung: Filterstrukturen ## +% ## -------------------------------------------------------------------------- ## +% ## Benoetigte(s) m-File(s): lcascade.m ## +% ################################################################################## + +% ##### Teilaufgabe b: Filterkoeffizienten berechnen ##### +[b a] = ellip(7,0.1,40,0.4); % Filter entwerfen +[B A] = lcascade(b,a); % 3-te kanonische Form +[Hg W] = freqz(b,a,512); % Uebertragungsfunktion des Gesamtsystems +[m n] = size(B); % Produkt der Uebertragungsfunktionen der Teilsysteme +Hi = zeros(512,m); + +for k=1:m; + Hi(:,k) = freqz(B(k,:),A(k,:),512); +end; +H3k = Hi(:,1) .* Hi(:,2) .* Hi(:,3) .* Hi(:,4); + +figure; +plot(W/pi,20*log10(abs(H3k)),'--b'); hold on; +plot(W/pi,20*log10(abs(Hg)),'-.b'); hold off; +grid; xlabel('Omega/pi (Nyquist = 1)'); ylabel('|H(exp(j*Omega))| in dB'); +title('-- Gesamtsystem, -. Produkt der Teilsysteme'); +% ##### EOF ##### diff --git a/common/l13.m b/common/l13.m new file mode 100755 index 0000000..4b70846 --- /dev/null +++ b/common/l13.m @@ -0,0 +1,47 @@ +% ################################################################################## +% ## Loesung: Entwurf rekursiver Filter ## +% ################################################################################## + +fD = 1000; %Durchlassfrequenz (Hz) +fS = 1400; %Sperrfrequenz (Hz) +Rp = 0.5; %Durchlaßripple (dB) +Rs = 30; %Sperrdaempfung (dB) +fA = 8000; %Abtastfrequenz (Hz) +Wp = 2*fD/fA; %Umrechnung in normierte Grenzfrequenzen +Ws = 2*fS/fA; + +for typ=1:4 + if typ == 1; + [n,Wn] = buttord(Wp,Ws,Rp,Rs); %Butterworth + [b,a] = butter(n,Wn); + FTyp = 'Butterworth'; + elseif typ == 2 ; + [n,Wn] = cheb1ord(Wp,Ws,Rp,Rs); %Tschebyscheff Typ I + [b,a] = cheby1(n,Rp,Wn); + FTyp = 'Tschebyscheff Typ I'; + elseif typ == 3 ; + [n,Wn] = cheb2ord(Wp,Ws,Rp,Rs); %Tschebyscheff Typ II + [b,a] = cheby2(n,Rs,Wn); + FTyp = 'Tschebyscheff Typ II'; + elseif typ == 4; + [n,Wn] = ellipord(Wp,Ws,Rp,Rs); %Cauer + [b,a] = ellip(n,Rp,Rs,Wn); + FTyp = 'Cauer'; + else; + error('Irgendwas stimmt nicht...'); + end; + + [H,W]=freqz(b,a,512); + + TitleString = sprintf('Filtertyp: %s Ordnung: %d', FTyp, n); + figure; plot(W/pi,abs(H)); grid; title(TitleString); + xlabel('Omega/pi');ylabel('|H(exp(j*Omega))|'); axis([0 1 0 1.2]); + xline=[Wp Wp];yline=[0 10.^(-Rp/20)]; line(xline,yline,'linestyle','--'); + xline=[0 Ws]; yline=[1 1]; line(xline,yline,'linestyle','--'); + xline=[0 Wp]; yline=[10.^(-Rp/20) 10.^(-Rp/20)];line(xline,yline,'linestyle','--'); + xline=[Ws Ws];yline=[10.^(-Rs/20) 1];line(xline,yline,'linestyle','--'); + xline=[Ws 1]; yline=[10.^(-Rs/20) 10.^(-Rs/20)];line(xline,yline,'linestyle','--'); + figure; zplane(b,a); title(TitleString); xlabel('Realteil'); + ylabel('Imaginaerteil'); +end; +% ##### EOF ##### diff --git a/common/l14.m b/common/l14.m new file mode 100755 index 0000000..53de299 --- /dev/null +++ b/common/l14.m @@ -0,0 +1,101 @@ +% ################################################################################## +% ## Loesung: Quantisierung der Filterkoeffizienten ## +% ## -------------------------------------------------------------------------- ## +% ## Benoetigte(s) m-File(s): lfxquant.m, lcoefrnd.m, lcascade.m ## +% ################################################################################## + +% ##### Teilaufgabe a: Cauer Tiefpass ##### +[b a] = ellip(7, 0.1, 40, 0.4); +[H W] = freqz(b,a,512); + +figure; plot(W/pi,20*log10(abs(H))); grid; xlabel('Omega/pi (Nyquist = 1)'); +ylabel('|H(exp(j*Omega))| in dB'); title('Cauer-TP'); + +% ##### Teilaufgabe b: Filterkoeffizienten quantisieren ##### +[b a] = ellip(7, 0.1, 40, 0.4); +[H W] = freqz(b,a,512); +[bq nb] = lcoefrnd(b,8); +[aq na] = lcoefrnd(a,8); +Hq = nb/na*freqz(bq,aq,512); + +figure; plot(W/pi,20*log10(abs(H))); hold on; plot(W/pi,20*log10(abs(Hq)),'-.'); +hold off; grid; title('- unquantisierte Koeff., -. 8 Bit quantisierte Koeff.'); +xlabel('Omega/pi (Nyquist = 1)'); ylabel('|H(exp(j*Omega))| in dB'); + +% ##### Teilaufgabe c: minimale Wortlaenge bestimmen ##### +[b a] = ellip(7,0.1,40,0.4); % Filter entwerfen +[H W] = freqz(b,a,512); % unquantisierte Uebertragungsfunktion +v_laenge=[16]; + +for l=v_laenge + + [bq,nb] = lcoefrnd(b,l); % Koeffizienten mit l Bits quantisieren + [aq,na] = lcoefrnd(a,l); + Hq = nb/na*freqz(bq,aq,512); % quantisierte Uebertragungsfunktion + figure; % Durchlassbereich darstellen + plot(W(1:208)/pi,20*log10(abs(H(1:208)))); axis([0 0.5 -1.4 0.3]); grid; hold on; + plot(W(1:208)/pi,-0.15*ones(size(Hq(1:208))),'w' ); %Toleranzgrenzen darstellen + plot(W(1:208)/pi,+0.15*ones(size(Hq(1:208))),'w' ); + plot(W(1:208)/pi,20*log10(abs(Hq(1:208))),'--'); hold off; + xlabel('Omega/pi (Nyquist = 1)');ylabel('|H(exp(j*Omega))| in dB'); + string = sprintf('Durch.: - unquantisiert, -- quantisiert mit %d Bits', l); + title(string); + figure; % Sperrbereich darstellen + plot(W(215:512)/pi,20*log10(abs(H(215:512)))); hold on; + plot(W(215:512)/pi,20*log10(abs(Hq(215:512))),'--'); + plot(W(215:512)/pi,-38*ones(size(Hq(215:512))),'w' ); grid; + hold off; xlabel('Omega/pi (Nyquist = 1)'); ylabel('|H(exp(j*Omega))| in dB'); + string = sprintf('Sperr.: - unquantisiert, -- quantisiert mit %d Bits', l); + title(string); + +end; + +% ##### Teilaufgabe d: Filterkoeffizienten quantisieren ##### + +[b a] = ellip(7,0.1,40,0.4); % Filter entwerfen +[B A] = lcascade(b,a); % unquantisierte Uebertragungsfunktion in + % 3-ter kanonischer Form berechnen +[m n] = size(B); % Produkt der Uebertragungsfunkt. der Teilsysteme +Hi = zeros(512,m); + +for k=1:m; + Hi(:,k) = freqz(B(k,:),A(k,:),512); +end; +H3k = Hi(:,1) .* Hi(:,2) .* Hi(:,3) .* Hi(:,4); +Bq = zeros(size(B)); % Vektoren initialisieren +Aq = Bq; +nb = zeros(m,1); +na = nb; +v_laenge=[12]; + +for l=v_laenge + + for k = 1:m; % Koeffizienten der Teilsysteme mit l Bits quantisieren + [Bq(k,:),nb(k)] = lcoefrnd(B(k,:),l); + [Aq(k,:),na(k)] = lcoefrnd(A(k,:),l); + end; + Hq = zeros(512,m); % quantisierte Gesamtuebertragungsfunktion + + for k=1:m + Hq(:,k) = nb(k)/na(k)*freqz(Bq(k,:),Aq(k,:),512); + end; + Hqg = Hq(:,1) .* Hq(:,2) .* Hq(:,3) .* Hq(:,4); + + figure; % Durchlassbereich darstellen + plot(W(1:208)/pi,20*log10(abs(H3k(1:208))));axis([0 0.5 -1.4 0.3]); grid; hold on; + plot(W(1:208)/pi,-0.15*ones(size(Hqg(1:208))),'w' ); %Toleranzgrenzen darstellen + plot(W(1:208)/pi,+0.15*ones(size(Hqg(1:208))),'w' ); + plot(W(1:208)/pi,20*log10(abs(Hqg(1:208))),'--'); hold off; + xlabel('Omega/pi (Nyquist = 1)'); ylabel('|H(exp(j*Omega))| in dB'); + string = sprintf('Durch.: - unquantisiert, -- quantisiert mit %d Bits', l); + title(string); + figure; % Sperrbereich darstellen + plot(W(215:512)/pi,20*log10(abs(H3k(215:512)))); hold on; + plot(W(215:512)/pi,20*log10(abs(Hqg(215:512))),'--'); + plot(W(215:512)/pi,-38*ones(size(Hqg(215:512))),'w' ); grid; hold off; + xlabel('Omega/pi (Nyquist = 1)'); ylabel('|H(exp(j*Omega))| in dB'); + string = sprintf('Sperr.: - unquantisiert, -- quantisiert mit %d Bits', l); + title(string); + +end; +% ##### EOF ##### diff --git a/common/l15.m b/common/l15.m new file mode 100755 index 0000000..7f866b0 --- /dev/null +++ b/common/l15.m @@ -0,0 +1,111 @@ +% ################################################################################## +% ## Loesung: Fensterfunktionen ## +% ## -------------------------------------------------------------------------- ## +% ## Benoetigte(s) m-File(s): ldtft.m ## +% ################################################################################## + + +% ##### Teilaufgabe a: Fenstertyp berechnen ##### +window_length=[21 42]; + +for n=window_length %%%Rechteckfenster%%% + n_find = find(n==window_length); + tltstr = 'Rechteck-Fenster'; + rw = boxcar(n); + [Rw W] = ldtft(rw,512); + + i=find(abs(Rw)==0); % Vermeidung von Werten==0 + Rw(i)=1e-100*ones(size(i)); + if n_find==1 + figure;plot(W/pi,20*log10(abs(Rw/max(Rw)))); + axis([-1 1 -80 0]);grid;hold on; + end; + plot(W/pi, 20*log10(abs(Rw/max(Rw))),'--'); xlabel('Omega/pi'); + ylabel('|F(exp(j*Omega))/F(exp(j*0))| in dB');title([tltstr,' m: - 21 -- 42']); +end; +hold off; +figure; stem(0:n-1,rw);grid; xlabel('k'); ylabel('f(k)'); title(tltstr); + +for n=window_length %%%von Hann%%% + n_find = find(n==window_length); + tltstr = 'von Hann-Fenster'; + rh = hanning(n); + [Rh W] = ldtft(rh,512); + + i=find(abs(Rh)==0); % Vermeidung von Werten==0 + Rh(i)=1e-100*ones(size(i)); + if n_find==1 + figure;plot(W/pi,20*log10(abs(Rh/max(Rh)))); axis([-1 1 -80 0]);grid;hold on; + end; + plot(W/pi,20*log10(abs(Rh/max(Rh))),'--'); xlabel('Omega/pi'); + ylabel('|F(exp(j*Omega))/F(exp(j*0))| in dB');title([tltstr,' m: - 21 -- 42']); +end; +hold off; +figure; stem(0:n-1,rh);grid; xlabel('k'); ylabel('f(k)'); title(tltstr); + +for n=window_length %%%Hamming%%% + n_find = find(n==window_length); + tltstr = 'Hamming-Fenster'; + rhm = hamming(n); + [Rhm W] = ldtft(rhm,512); + + i=find(abs(Rhm)==0); % Vermeidung von Werten==0 + Rhm(i)=1e-100*ones(size(i)); + if n_find==1 + figure;plot(W/pi,20*log10(abs(Rhm/max(Rhm))));axis([-1 1 -80 0]);grid;hold on; + end; + plot(W/pi,20*log10(abs(Rhm/max(Rhm))),'--'); xlabel('Omega/pi'); + ylabel('|F(exp(j*Omega))/F(exp(j*0))| in dB');title([tltstr,' m: - 21 -- 42']); +end; +hold off; +figure; stem(0:n-1,rhm);grid; xlabel('k'); ylabel('f(k)'); title(tltstr); + +for n=window_length %%%Blackman%%% + n_find = find(n==window_length); + tltstr = 'Blackman-Fenster'; + rhb = blackman(n); + [Rhb W] = ldtft(rhb,512); + + i=find(abs(Rhb)==0); % Vermeidung von Werten==0 + Rhb(i)=1e-100*ones(size(i)); + if n_find==1 + figure;plot(W/pi,20*log10(abs(Rhb/max(Rhb))));axis([-1 1 -80 0]);grid;hold on; + end; + plot(W/pi,20*log10(abs(Rhb/max(Rhb))),'--'); xlabel('Omega/pi'); + ylabel('|F(exp(j*Omega))/F(exp(j*0))| in dB');title([tltstr,' m: - 21 -- 42']); +end; +hold off; +figure; stem(0:n-1,rhb);grid; xlabel('k'); ylabel('f(k)'); title(tltstr); + +for n=window_length %%%Bartlett%%% + n_find = find(n==window_length); + tltstr = 'Bartlett-Fenster'; + rhbt = bartlett(n); + [Rhbt W] = ldtft(rhbt,512); + + i=find(abs(Rhbt)==0); % Vermeidung von Werten==0 + Rhbt(i)=1e-100*ones(size(i)); + if n_find==1 + figure;plot(W/pi,20*log10(abs(Rhbt/max(Rhbt))));axis([-1 1 -80 0]);grid; + hold on; + end; + plot(W/pi,20*log10(abs(Rhbt/max(Rhbt))),'--'); xlabel('Omega/pi'); + ylabel('|F(exp(j*Omega))/F(exp(j*0))| in dB');title([tltstr,' m: - 21 -- 42']); +end; +hold off; +figure; stem(0:n-1,rhbt);grid; xlabel('k'); ylabel('f(k)'); title(tltstr); + +% ##### Teilaufgabe b: Kaiser-Fenster berechnen ##### +w0 = kaiser(21,0); [W0 W] = ldtft(w0,512); +w3 = kaiser(21,3); W3 = ldtft(w3,512); +w6 = kaiser(21,6); W6 = ldtft(w6,512); +k=0:20; + +figure; plot(W/pi,20*log10(abs(W0/max(W0)))); xlabel('Omega/pi');axis([-1 1 -80 0]); +ylabel('|F(exp(j*Omega))/F(exp(j*0))| in dB'); hold on; +plot(W/pi,20*log10(abs(W3/max(W3))),'--');plot(W/pi,20*log10(abs(W6/max(W6))),'-.'); +hold off; title('Kaiser-Fenster: beta: - 0, -- 3, -. 6'); grid; +figure; plot(k,w0); xlabel('k');ylabel('f(k)');hold on; plot(k,w3,'--'); +plot(k,w6,'-.'); axis([0 20 0 1.2]); hold off; +title('Kaiser-Fenster: beta: - 0, -- 3, -. 6 '); grid +% ##### EOF ##### diff --git a/common/l16.m b/common/l16.m new file mode 100755 index 0000000..4bf5565 --- /dev/null +++ b/common/l16.m @@ -0,0 +1,35 @@ +% ################################################################################## +% ## Loesung: Entwurf nichtrekursiver Filter mittels Fensterung ## +% ################################################################################## + +% ##### Teilaufgabe a: Tiefpass entwerfen ##### +bR = fir1(50,0.4,boxcar(51)); [HR W] = freqz(bR,1,512); +bH = fir1(50,0.4,hanning(51)); [HH W] = freqz(bH,1,512); +bHm = fir1(50,0.4,hamming(51)); [HHM W] = freqz(bHm,1,512); +bB = fir1(50,0.4,blackman(51)); [HB W] = freqz(bB,512); +bK = fir1(50,0.4,kaiser(51,8)); [HK W] = freqz(bK,512); + +plot(W/pi,20*log10(abs(HR/max(HR)))); hold on; +plot(W/pi,20*log10(abs(HH/max(HH))),'--'); +plot(W/pi,20*log10(abs(HHM/max(HHM))),'-.'); hold off; +xlabel('Omega/pi'); ylabel('|H(exp(j*Omega))/H(1)| in dB'); grid +title(['Tiefpass: - Rechteck -- von Hann -. Hamming']); axis([0 1 -150 0]); figure; +plot(W/pi,20*log10(abs(HB/max(HB))),'-'); hold on; +plot(W/pi,20*log10(abs(HK/max(HK))),'--'); hold off; +xlabel('Omega/pi'); ylabel('|H(exp(j*Omega))/H(1)| in dB'); grid +title(['Tiefpass: - Blackman -- Kaiser mit beta=8']); axis([0 1 -150 0]); + +% ##### Teilaufgabe b: Hochpass entwerfen ##### +b = fir1(33,0.4,'high'); [H W] = freqz(b,1,512); + +figure; +plot(W/pi,20*log10(abs(H/max(H)))); grid; +xlabel('Omega/pi'); ylabel('|H(exp(j*Omega))/H(-1)| in dB'); title('Hochpass'); + +% ##### Teilaufgabe c: Gruppenlaufzeit bestimmen ##### +figure; +b = fir1(50,0.4); +grpdelay(b,1,256); +title('Gruppenlaufzeit in Abtastwerten'); ylabel(''); xlabel('Omega/pi'); +axis([0 1 24 26]); +% ##### EOF ##### diff --git a/common/l17.m b/common/l17.m new file mode 100755 index 0000000..9097c91 --- /dev/null +++ b/common/l17.m @@ -0,0 +1,51 @@ +% ################################################################################## +% ## Loesung: Remez-Entwurf (Tschebyscheff-Approximation) ## +% ################################################################################## + +% ##### Teilaufgabe a: Filter entwerfen ##### +f = [0 0.3 0.5 1.0]; +g = [1 1 0 0]; +v_m=[9 10]; % Vektor mit den Filterordnungen + +for m=v_m + b = remez(m,f,g); [H W] = freqz(b,1,512); + string = sprintf('Gewuenschte/Aktuelle Antwort - Filterordnung: %d', m); + figure;plot(f,g,W/pi,abs(H)); + line([0 0.3],[ 2-max(abs(H)) 2-max(abs(H))],'linestyle','--'); + line([0 0.5],[max(abs(H)) max(abs(H))],'linestyle','--'); + line([0.5 1],[max(abs(H(400:length(H)))) max(abs(H(400:length(H))))],... + 'linestyle','--'); + line([0.3 0.3],[ 2-max(abs(H)) 0],'linestyle','--'); + line([0.5 0.5],[ max(abs(H)) max(abs(H(400:length(H))))],'linestyle','--'); + grid; xlabel('Omega/pi'); ylabel('|H(exp(j*Omega))|'); title(string); +end; + +% ##### Teilaufgabe b: Frequenzgang wichten ##### +f = [0 0.3 0.5 1.0]; +g = [1 1 0 0]; +w = [ 0.5 2 ]; +v_m=[9 10]; % Vektor mit den Filterordnungen + +for m=v_m + b = remez(m,f,g,w); [H W] = freqz(b,1,512); + string = sprintf('Gewuenschte/Aktuelle Antwort - Filterordnung: %d', m); + figure; plot(f,g,W/pi,abs(H)); + line([0 0.3],[ 2-max(abs(H)) 2-max(abs(H))],'linestyle','--'); + line([0 0.5],[max(abs(H)) max(abs(H))],'linestyle','--'); + line([0.5 1],[max(abs(H(400:length(H)))) max(abs(H(400:length(H))))],... + 'linestyle','--'); + line([0.3 0.3],[ 2-max(abs(H)) 0],'linestyle','--'); + line([0.5 0.5],[ max(abs(H)) max(abs(H(400:length(H))))],'linestyle','--'); + grid; xlabel('Omega/pi'); ylabel('|H(exp(j*Omega))|'); title(string); +end; + +% ##### Teilaufgabe c: Tiefpass-Bandpass Kombination ##### +f = [0 0.2 0.26 0.44 0.5 0.7 0.76 1]; % Vektor der Eckfrequenzen +g = [1 1 0 0 1 1 0 0]; % Vektor der Frequenzgaenge +w = [ 1 4 1 4]; % Vektor mit den Gewichten +m = 128; % Filterordnung +b = remez(m,f,g,w); [H W] = freqz(b,1,512); + +figure; plot(W/pi,20*log10(abs(H)) );grid; xlabel('Omega/pi'); +ylabel('|H(exp(j*Omega))| in dB'); title('TP-BP-Kombination'); +% ##### EOF ##### diff --git a/common/l18.m b/common/l18.m new file mode 100755 index 0000000..60c02ba --- /dev/null +++ b/common/l18.m @@ -0,0 +1,31 @@ +% ################################################################################## +% ## Loesung: Entwurf eines Differenzieres (Remez-Verfahren) ## +% ################################################################################## + +% ##### Teilaufgabe a: Differenzierer entwerfen ##### +% ### 1: +f = [0 1]; m = [0 1]; b1 = remez(22,f,m); +f = [0 0.9]; m = [0 1]; b2 = remez(21,f,m); + +figure; [H1 W] = freqz(b1,1,512); plot(W/pi,abs(H1) ); grid; xlabel('Omega/pi'); +ylabel('|H(exp(j*Omega))|'); title('Filterordnung: 22, Typ I'); +figure; [H2 W] = freqz(b2,1,512); plot(W/pi,abs(H2) ); grid; xlabel('Omega/pi'); +ylabel('|H(exp(j*Omega))|'); title('Filterordnung: 21, Typ III'); +figure; stem([0:22], b1); axis([0 22 -0.3 0.5]); grid; xlabel('k'); +ylabel('Impulsantwort'); title('Filterordnung: 22, Typ I'); +figure; stem([0:21], b2); axis([0 21 -0.3 0.5]); grid; xlabel('k'); +ylabel('Impulsantwort'); title('Filterordnung: 21, Typ III'); + +% ##### Teilaufgabe b: Uebertragungsfunktion vergleichen ##### +f = [0 1]; m = [0 1]; b2 = remez(21,f,m,'d'); +f = [0 0.9]; m = [0 1]; b1 = remez(22,f,m,'d'); + +figure; [H1 W] = freqz(b1,1,512); plot(W/pi,abs(H1) ); grid; xlabel('Omega/pi'); +ylabel('|H_D(exp(j*Omega))|'); title('Filterordnung: 22, Typ II'); +figure; [H2 W] = freqz(b2,1,512); plot(W/pi,abs(H2)); grid; xlabel('Omega/pi'); +ylabel('|H_D(exp(j*Omega))|'); title('Filterordnung: 21, Typ IV'); +figure; stem([0:22], b1); axis([0 22 -0.4 0.5]); grid; xlabel('k'); +ylabel('Impulsantwort'); title('Filterordnung: 22, Typ II'); +figure;stem([0:21], b2); axis([0 21 -0.4 0.5]); grid; xlabel('k'); +ylabel('Impulsantwort'); title('Filterordnung: 21, Typ IV'); +% ##### EOF ##### diff --git a/common/l19.m b/common/l19.m new file mode 100755 index 0000000..f62dda8 --- /dev/null +++ b/common/l19.m @@ -0,0 +1,40 @@ +% ################################################################################## +% ## Loesung: Entwurf eines Hilbert-Transformators (Remez-Verfahren) ## +% ################################################################################## + +% ##### Teilaufgabe a: Hilbert-Transformator entwerfen ##### +b1 = remez(30,[0.05 0.95],[1 1],'h'); +b2 = remez(31,[0.05 1],[1 1],'h'); + +figure; [H1 W] = freqz(b1,1,512);plot(W/pi,abs(H1)); grid; xlabel('Omega/pi'); +ylabel('|H_H(exp(j*Omega))|'); title('Filterordnung: 30 Typ II'); +figure; [H2 W] = freqz(b2,1,512);plot(W/pi,abs(H2)); grid; xlabel('Omega/pi'); +ylabel('|H_H(exp(j*Omega))|'); title('Filterordnung: 31 Typ IV'); +figure; stem([0:30], b1); axis([0 30 -0.8 0.8]); grid; xlabel('k'); +ylabel('Impulsantwort'); title('Filterordnung: 30 Typ II'); +figure; stem([0:31], b2); axis([0 31 -0.8 0.8]); grid; xlabel('k'); +ylabel('Impulsantwort'); title('Filterordnung: 31 Typ IV'); + +% ##### Teilaufgabe b: Ausgangssignal filtern ##### +b = remez(30,[0.05 0.95],[1 1],'h'); +Fs = 1000; +t = (0:1/Fs:2)'; +x = sin(2*pi*50*t); +xh = filter(b,1,x); + +figure, zplane(b); xlabel('Realteil'); ylabel('Imaginaerteil'); +title('Pol-Nullstellen-Diagramm'); +xd = [zeros(15,1); x(1:length(x)-15)]; % delay 15 samples +figure; plot(t(1:50),xd(1:50),t(1:50),xh(1:50),'--'); grid; xlabel('Zeit (sek)'); +ylabel('Amplitude'); title(' - Eingangssig. -- Ausgangssig.'); + +% ##### Teilaufgabe c: Hilbert-Tranformierte berechnen ##### +Fs = 1000; +t = (0:1/Fs:2)'; +x = sin(2*pi*50*t); +y = hilbert(x); + +figure; plot(t(1:50), real(y(1:50))); grid; hold on; +plot(t(1:50), imag(y(1:50)),'--'); xlabel('Zeit (sek)'); ylabel('Amplitude'); +title('-- Hilbert-Transform. - Original'); +% ##### EOF ##### diff --git a/common/l2.m b/common/l2.m new file mode 100755 index 0000000..a3aa3b3 --- /dev/null +++ b/common/l2.m @@ -0,0 +1,48 @@ +% ################################################################################## +% ## Loesung: Trigonometrische Folge ## +% ## -------------------------------------------------------------------------- ## +% ## Benoetigte(s) m-File(s): ldsin.m ## +% ################################################################################## + +% ##### Teilaufgabe a: Graphische Darstellung von Signalen ##### +% ### 1: +k = 0:25; +x1 = sin((pi/17)*k); + +figure; stem(k,x1);grid; xlabel('k'); ylabel('x1(k)'); + +% ### 2: +k = -15:25; +x2 = sin((pi/17)*k); + +figure; stem(k,x2); grid; xlabel('k'); ylabel('x2(k)'); +set(gca,'XLim',[-15 25]); set(gca,'XTick',[-15 -10 -5 0 5 10 15 20 25]); + +% ### 3: +k = -10:10; +x3 = sin(3*pi*k + (pi/2) ); + +figure; stem(k,x3); grid; xlabel('k'); ylabel('x3(k)'); + +% ### 4: +k = 0:50; +x4 = cos( pi/sqrt(23) * k ); + +figure; stem(k,x4); grid; xlabel('k'); ylabel('x4(k)'); + +% ### 5: +k = 0:20; +x5 = cos( k ); + +figure; stem(k,x5); grid; xlabel('k'); ylabel('x5(k)'); + +% ##### Teilaufgabe b: function m-file ##### +[y,k] = ldsin(2,pi/11,pi/2,-20,20); + +figure; stem(k,y); grid; xlabel('k'); ylabel('x(k)'); + +% ##### Teilaufgabe c: Modifizierung ##### +[y,k] = ldsin(2,pi/11,pi/2,-20,20); + +figure; stem(k,y); grid; title('Teilaufgabe c'); xlabel('k'); ylabel('Amplitude'); +% ##### EOF ##### diff --git a/common/l20.m b/common/l20.m new file mode 100755 index 0000000..d871759 --- /dev/null +++ b/common/l20.m @@ -0,0 +1,39 @@ +% ################################################################################## +% ## Loesung: Diskrete Fouriertransformation (DFT) ## +% ################################################################################## + +% ##### Teilaufgabe a: DFT von x(k) berechnen ##### +kk = 0:15; +x = kk + j*(7-kk); xr=real(x); xi=imag(x); +X = fft(x); Xr=real(X); Xi=imag(X); +figure; stem(kk,Xr); title('Re(X(n))'); xlabel('n'); grid; +figure; stem(kk,Xi); title('Im(X(n))'); xlabel('n'); grid; +input('....Fuer naechste Teilaufgabe: RETURN druecken'); + +% ##### Teilaufgabe b: DFT von Re{x(k)} und Im{x(k)} ##### +figure; stem(kk,real(fft(xr))); title('Re(FFT(Re(x(k))))'); xlabel('n'); grid; +figure; stem(kk,imag(fft(xr))); title('Im(FFT(Re(x(k))))'); xlabel('n'); grid; +figure; stem(kk,real(fft(j*xi))); title('Re(FFT(j*Im(x(k))))'); xlabel('n'); grid; +figure; stem(kk,imag(fft(j*xi))); title('Im(FFT(j*Im(x(k))))'); xlabel('n'); grid; +input('....Fuer naechste Teilaufgabe: RETURN druecken'); + +% ##### Teilaufgabe c: Zyklische Verschiebung ##### +xtmp = [x x]; +x1_04 = xtmp(kk+5); X1_04 = fft(x1_04); +x1_08 = xtmp(kk+9); X1_08 = fft(x1_08); +figure; stem(kk,Xr); title('Re(X(n))'); xlabel('n'); grid; +figure; stem(kk,Xi); title('Im(X(n))'); xlabel('n'); grid; +figure; stem(kk,real(X1_04)); title('Re(X1(n)), lambda=4'); xlabel('n'); grid; +figure; stem(kk,imag(X1_04)); title('Im(X1(n)), lambda=4'); xlabel('n'); grid; +figure; stem(kk,real(X1_08)); title('Re(X1(n)), lambda=8'); xlabel('n'); grid; +figure; stem(kk,imag(X1_08)); title('Im(X1(n)), lambda=8'); xlabel('n'); grid; +input('....Fuer naechste Teilaufgabe: RETURN druecken'); + +% ##### Teilaufgabe d: Multiplikation mit exp(j*...) ##### +x2 = x.*exp(-j*2*pi*3*kk/16); +X2 = fft(x2); +figure; stem(kk,Xr); title('Re(X(n))'); xlabel('n'); grid; +figure; stem(kk,Xi); title('Im(X(n))'); xlabel('n'); grid; +figure; stem(kk,real(X2)); title('Re(X2(n))'); xlabel('n'); grid; +figure; stem(kk,imag(X2)); title('Im(X2(n))'); xlabel('n'); grid; +% ##### EOF ##### diff --git a/common/l21.m b/common/l21.m new file mode 100755 index 0000000..79ea5a3 --- /dev/null +++ b/common/l21.m @@ -0,0 +1,13 @@ +% ################################################################################## +% ## Loesung: DFT: Interpolation ## +% ################################################################################## + +N = 64; +x = 1-cos(2*pi/N*(0:N-1)); X = fft(x); +x0 = [x zeros(1,N)]; X0 = fft(x0); + +figure; stem(0:N-1,abs(X));axis([0 N-1 0 max(X)]);title('abs(X(n))'); +xlabel('n'); +figure; stem(0:2*N-1,abs(X0)); axis([0 2*N-1 0 max(X)]);title('abs(X0(n))'); +xlabel('n'); +% ##### EOF ##### diff --git a/common/l22.m b/common/l22.m new file mode 100755 index 0000000..81ad2e5 --- /dev/null +++ b/common/l22.m @@ -0,0 +1,16 @@ +% ################################################################################## +% ## Loesung: Alternative Berechnung der IFFT ## +% ################################################################################## + +N = 64; +x = 1 - cos(2*pi/N*(0:N-1)); +x0 = [x zeros(1,N)]; + +X0 = fft(x0); % FFT +X1 = imag(X0) + j*real(X0); % Re() und Im() vertauschen +x1 = fft(X1); % Nochmalige FFT +x2 = imag(x1) + j*real(x1); % Re() und Im() vertauschen + +figure; bar(0:2*N-.5,x0); title('x0(k)'); xlabel('k'); axis([-.5 2*N-1 0 2.1]); +figure; bar(0:2*N-.5,x2); title('x2(k)'); xlabel('k'); axis([-.5 2*N-1 0 269]); +% ##### EOF ##### diff --git a/common/l23.m b/common/l23.m new file mode 100755 index 0000000..ef1ae2c --- /dev/null +++ b/common/l23.m @@ -0,0 +1,12 @@ +% ################################################################################## +% ## Loesung: DFT reeller Datenfolgen ## +% ## -------------------------------------------------------------------------- ## +% ## Benoetigte(s) m-File(s): lfft.m ## +% ################################################################################## + +N = 2^10; +x = randn(1,N); +X1 = fft(x); +X2 = lfft(x); +differenz = mean(abs(X1-X2).^2) +% ##### EOF ##### diff --git a/common/l24.m b/common/l24.m new file mode 100755 index 0000000..2198eca --- /dev/null +++ b/common/l24.m @@ -0,0 +1,19 @@ +% ################################################################################## +% ## Loesung: Spektraltransformation reeller Bandpass-Signale ## +% ## -------------------------------------------------------------------------- ## +% ## Benoetigte(s) m-File(s): lband.m ## +% ################################################################################## + +NFFT = 2^10; +b1 = lband(1000,0.2,0.02); fT = 0:1/NFFT:1-1/NFFT; + % Originalsignal +b2 = b1; lb = length(b2); aB2_hn=abs(fft(b2.*hanning(lb),NFFT)); +figure; plot(fT, aB2_hn); axis([0 1 0 1.1*max(aB2_hn)]); xlabel('Omega/2pi = f T'); +ylabel('|B(exp(j*Omega))|'); + +for k=[9 10 11 12] %Unterabtastung um Faktor k + b2 = b1(1:k:1000); lb = length(b2); aB2_hn=abs(fft(b2.*hanning(lb),NFFT)); + figure; plot(fT,aB2_hn); axis([0 1 0 1.1*max(aB2_hn)]); xlabel('Omega/2pi = f T'); + ylabel('|B_M(exp(j*Omega))|'); title(sprintf('Unterabtastung um Faktor M=%d',k)); +end; +% ##### EOF ##### diff --git a/common/l25.m b/common/l25.m new file mode 100755 index 0000000..fcfc0df --- /dev/null +++ b/common/l25.m @@ -0,0 +1,31 @@ +% ################################################################################## +% ## Loesung: DFT: Abbruchfehler, Unterabtastung ## +% ################################################################################## + +kk08 = 0:7; +kk64 = 0:63; +x1_08 = sin(kk08*pi/8); x1_64 = sin(kk64*pi/8); +x2_08 = sin(kk08*pi/2); x2_64 = sin(kk64*pi/2); +x3_08 = sin(kk08*pi*3/2); x3_64 = sin(kk64*pi*3/2); + +figure; kk100=0:8/100:8; plot(kk100,sin(kk100*pi/8),':'); hold on; +stem(kk08,x1_08); axis([0 8 -1.2 1.2]); title('x1(k)=sin(k*pi/8), k=0...7'); +figure;stem(kk08,abs(fft(x1_08))); axis([0 7 -1 6]); title('abs(X1(n)),n=0...7'); +figure;stem(kk64,x1_64);axis([0 64 -1.2 1.2]);title('x1(k)=sin(k*pi/8),k=0...63'); +figure;stem(kk64,abs(fft(x1_64)));axis([0 64 -1 40]);title('abs(X1(n)),n=0...63'); +input('....Zum Fortfahren: RETURN druecken'); + +figure; kk100=0:8/100:8; plot(kk100,sin(kk100*pi/2),':'); hold on; +stem(kk08,x2_08); axis([0 8 -1.2 1.2]); title('x2(k)=sin(k*pi/2), k=0...7'); +figure;stem(kk08,abs(fft(x2_08))); axis([0 7 -1 6]); title('abs(X2(n)),n=0...7'); +figure;stem(kk64,x2_64);axis([0 64 -1.2 1.2]);title('x2(k)=sin(k*pi/2),k=0...63'); +figure;stem(kk64,abs(fft(x2_64)));axis([0 64 -1 40]);title('abs(X2(n))),n=0...63'); +input('....Zum Fortfahren: RETURN druecken'); + +figure; kk100=0:8/100:8; plot(kk100,sin(kk100*pi*3/2),':'); hold on; +stem(kk08,x3_08); axis([0 8 -1.2 1.2]); title('x3(k)=sin(k*pi*3/2), k=0...7'); +figure;stem(kk08,abs(fft(x3_08))); axis([0 8 -1 6]); title('abs(X3(n)), n=0...7'); +axis([0 7 -1 6]); +figure;stem(kk64,x3_64);axis([0 64 -1.2 1.2]);title('x3(k)=sin(k*pi*3/2),k=0...63'); +figure;stem(kk64,abs(fft(x3_64)));axis([0 64 -1 40]);title('abs(X3(n)),n=0...63'); +% ##### EOF ##### diff --git a/common/l26.m b/common/l26.m new file mode 100755 index 0000000..f97b5d9 --- /dev/null +++ b/common/l26.m @@ -0,0 +1,69 @@ +% ################################################################################## +% ## Loesung: Vergleich versch. Fensterfunktionen ## +% ################################################################################## +N = 63; NFFT = 2^9; +f_re = boxcar(N); aF_re = abs(fft(f_re,NFFT)/sum(f_re)); +f_hn = hanning(N); aF_hn = abs(fft(f_hn,NFFT)/sum(f_hn)); +f_hm = hamming(N); aF_hm = abs(fft(f_hm,NFFT)/sum(f_hm)); +f_bl = blackman(N); aF_bl = abs(fft(f_bl,NFFT)/sum(f_bl)); +a_d1=48.6; f_d1=chebwin(N,a_d1); aF_d1 = abs(fft(f_d1,NFFT)/sum(f_d1)); +a_d2=76.1; f_d2=chebwin(N,a_d2); aF_d2 = abs(fft(f_d2,NFFT)/sum(f_d2)); +kk = 0:N-1; fT = 0:1/NFFT:1-1/NFFT; +% Berechnung von a_{min} (fsT:=Omega_s/(2*pi); Sperrber.: n_s+1<=n<=NFFT+1-n_s) +fsT=1/N; n_s=round(NFFT*fsT); a_re=20*log10(max(aF_re(n_s+1:NFFT+1-n_s))); +fsT=2/N; n_s=round(NFFT*fsT); a_hn=20*log10(max(aF_hn(n_s+1:NFFT+1-n_s))); +fsT=2/N; n_s=round(NFFT*fsT); a_hm=20*log10(max(aF_hm(n_s+1:NFFT+1-n_s))); +fsT=3/N; n_s=round(NFFT*fsT); a_bl=20*log10(max(aF_bl(n_s+1:NFFT+1-n_s))); + +figure; % RECHTECK-Fenster +bar(kk,f_re); axis([-.5 N-.5 0 1.1]); title('Rechteck-Fenster'); ylabel('f(k)'); +xlabel('k'); +figure; plot(fT,20*log10(aF_re)); hold on; axis([0 1 -150 10]); +text(.3,a_re+13,sprintf('a_{min} = %3.1f dB',-a_re)); +plot([0 1.9],[a_re a_re],'--'); title('|Normierte Spektralfunktion|'); +ylabel('|F(exp(j*Om.))/F(exp(j*0))| in dB'); xlabel('Omega/2pi = f T'); +input('....Zum Fortfahren: RETURN druecken'); + +figure; % HANNING-Fenster +bar(kk,f_hn); axis([-.5 N-.5 0 1.1]); title('Hanning-Fenster'); ylabel('f(k)'); +xlabel('k'); +figure; plot(fT,20*log10(aF_hn)); hold on; axis([0 1 -150 10]); +text(.3,a_hn+15,sprintf('a_{min} = %3.1f dB', -a_hn)); +plot([0 1.9],[a_hn a_hn],'--'); title('|Normierte Spektralfunktion|'); +ylabel('|F(exp(j*Om.))/F(exp(j*0))| in dB'); xlabel('Omega/2pi = f T'); +input('....Zum Fortfahren: RETURN druecken'); + +figure; % HAMMING-Fenster +bar(kk,f_hm); axis([-.5 N-.5 0 1.1]); title('Hamming-Fenster'); ylabel('f(k)'); +xlabel('k'); +figure; plot(fT,20*log10(aF_hm)); hold on; axis([0 1 -150 10]); +text(.3,a_hm+15,sprintf('a_{min} = %3.1f dB', -a_hm)); +plot([0 1.9],[a_hm a_hm],'--'); title('|Normierte Spektralfunktion|'); +ylabel('|F(exp(j*Om.))/F(exp(j*0))| in dB'); xlabel('Omega/2pi = f T'); +input('....Zum Fortfahren: RETURN druecken'); + +figure; % BLACKMAN-Fenster +bar(kk,f_bl); axis([-.5 N-.5 0 1.1]); title('Blackman-Fenster'); ylabel('f(k)'); +xlabel('k'); +figure; plot(fT,20*log10(aF_bl)); hold on; axis([0 1 -150 10]); +text(.3,a_bl+15,sprintf('a_{min} = %3.1f dB', -a_bl)); +plot([0 1.9],[a_bl a_bl],'--'); title('|Normierte Spektralfunktion|'); +ylabel('|F(exp(j*Om.))/F(exp(j*0))| in dB'); xlabel('Omega/2pi = f T'); +input('....Zum Fortfahren: RETURN druecken'); + +figure; % DOLPH-TSCHEBYSCHEFF-1 +bar(kk,f_d1); axis([-.5 N-.5 0 1.1]);title(sprintf('Dolph-Tscheby.,%3.1fdB',a_d1)); +ylabel('f(k)'); xlabel('k'); +figure; plot(fT,20*log10(aF_d1)); hold on; axis([0 1 -150 10]); +text(.3,-a_d1+15,sprintf('a_{min} = %3.1f dB',a_d1)); +title('|Normierte Spektralfunktion|'); ylabel('|F(exp(j*Om.))/F(exp(j*0))| in dB'); +xlabel('Omega/2pi = f T'); input('....Zum Fortfahren: RETURN druecken'); + +figure; % DOLPH-TSCHEBYSCHEFF-2 +bar(kk,f_d2); axis([-.5 N-.5 0 1.1]);title(sprintf('Dolph-Tscheby.,%3.1fdB',a_d2)); +ylabel('f(k)'); xlabel('k'); +figure; plot(fT,20*log10(aF_d2)); hold on; axis([0 1 -150 10]); +text(.3,-a_d2+15,sprintf('a_{min} = %3.1f dB',a_d2)); +title('|Normierte Spektralfunktion|'); ylabel('|F(exp(j*Om.))/F(exp(j*0))| in dB'); +xlabel('Omega/2pi = f T'); +% ##### EOF ##### diff --git a/common/l27.m b/common/l27.m new file mode 100755 index 0000000..f5d5be5 --- /dev/null +++ b/common/l27.m @@ -0,0 +1,33 @@ +% ################################################################################## +% ## Loesung: Seitenbanddaempfung des Blackman-Fenster ## +% ################################################################################## + +NFFT=2^10; +f_bl032= blackman(32); aF_bl032= abs(fft(f_bl032,NFFT)/sum(f_bl032)); +f_bl064= blackman(64); aF_bl064= abs(fft(f_bl064,NFFT)/sum(f_bl064)); +f_bl128= blackman(128); aF_bl128= abs(fft(f_bl128,NFFT)/sum(f_bl128)); + fT = 0:1/NFFT:1-1/NFFT; +% Berechnung der minimalen Sperrdaempfung a_{min}: +% Os := Omega_s/(2*pi); n_s: Sperrbereich-Indices: n_s+1<=n<=NFFT+1-n_s +Os=3/32; n_s=round(NFFT*Os); a032=20*log10(max(aF_bl032(n_s+1:NFFT+1-n_s))); +Os=3/64; n_s=round(NFFT*Os); a064=20*log10(max(aF_bl064(n_s+1:NFFT+1-n_s))); +Os=3/128; n_s=round(NFFT*Os); a128=20*log10(max(aF_bl128(n_s+1:NFFT+1-n_s))); + +figure; plot(fT, 20*log10(aF_bl032)); hold on; axis([0 1 -180 10]); +text(0.41,a032,sprintf('a_{min} = %3.1f dB',-a032)); +plot([0 0.39], [a032 a032], '--'); +title('|Norm. Spektralfunktion| des Blackmanfensters: N=32'); +ylabel('|F(exp(j*Om.))/F(exp(j*0))| in dB'); xlabel('Omega/2pi = f T'); + +figure; plot(fT, 20*log10(aF_bl064)); hold on; axis([0 1 -180 10]); +text(0.41,a064,sprintf('a_{min} = %3.1f dB',-a064)); +plot([0 0.39], [a064 a064], '--'); +title('|Norm. Spektralfunktion| des Blackmanfensters: N=64'); +ylabel('|F(exp(j*Om.))/F(exp(j*0))| in dB'); xlabel('Omega/2pi = f T'); + +figure; plot(fT, 20*log10(aF_bl128)); hold on; axis([0 1 -180 10]); +text(0.41,a128,sprintf('a_{min} = %3.1f dB',-a128)); +plot([0 0.39], [a128 a128], '--'); +title('|Norm. Spektralfunktion| des Blackmanfensters: N=128'); +ylabel('|F(exp(j*Om.))/F(exp(j*0))| in dB'); xlabel('Omega/2pi = f T'); +% ##### EOF ##### diff --git a/common/l28.m b/common/l28.m new file mode 100755 index 0000000..5ef9b88 --- /dev/null +++ b/common/l28.m @@ -0,0 +1,10 @@ +% ################################################################################## +% ## Loesung: Dolph-Tschebyscheff-Fenster ## +% ## -------------------------------------------------------------------------- ## +% ## Benoetigte(s) m-File(s): lcheby.m ## +% ################################################################################## +% +[f,F,a_min]=lcheby(32, 0.1*pi, 2*pi/256); +[f,F,a_min]=lcheby(64, 0.1*pi, 2*pi/256); +[f,F,a_min]=lcheby(32, 0.4*pi, 2*pi/256); +% ##### EOF ##### diff --git a/common/l29.m b/common/l29.m new file mode 100755 index 0000000..9cd4c66 --- /dev/null +++ b/common/l29.m @@ -0,0 +1,32 @@ +% ################################################################################## +% ## Loesung: Verminderung des Leckeffektes durch Fensterung ## +% ################################################################################## + +N = 64; +kk = (0:N-1); kk2 = (0:2*N-1); +x1 = sin(kk*pi/8); X1 = fft(x1); +x1Hm = x1.*hamming(N).'; X1Hm = fft(x1Hm); +x2 = sin(kk*pi/10); X2 = fft(x2); +x2Hm = x2.*hamming(N).'; X2Hm = fft(x2Hm); + +figure; stem(kk2, [x1 x1]); axis([0 2*N -1.1 1.1]); ylabel('x1((k))'); +title('x1(k)=sin(k*pi/8), k=0..63, zweifach periodisch fortgesetzt'); xlabel('k'); +figure; stem(kk, abs(X1)); axis([0 N 0 1.1*N/2]); ylabel('|X1(n)|'); +title('Betragsspektrum von x1(k)'); xlabel('n'); + +figure; stem(kk2, [x1Hm x1Hm]); axis([0 2*N -1.1 1.1]); ylabel('x1_Hm((k))'); +title('Hamming-gefenstertes x1(k), zweifach periodisch fortgesetzt'); xlabel('k'); +figure; stem(kk, abs(X1Hm)); axis([0 N 0 1.1*N/2]); ylabel('|X1_Hm(n)|'); +title('Betragsspektrum bei Hamming-Fensterung'); xlabel('n'); +input('......fuer 2.Datenfolge: RETURN druecken'); + +figure; stem(kk2, [x2 x2]); axis([0 2*N -1.1 1.1]); ylabel('x2((k))'); +title('x2(k)=sin(k*pi/10), k=0..63, zweifach periodisch fortgesetzt'); xlabel('k'); +figure; stem(kk, abs(X2)); axis([0 N 0 1.1*N/2]); ylabel('|X2(n)|'); +title('Betragsspektrum von x2(k)'); xlabel('n'); + +figure; stem(kk2, [x2Hm x2Hm]); axis([0 2*N -1.1 1.1]); ylabel('x2_Hm((k))'); +title('Hamming-gefenstertes x2(k), zweifach periodisch fortgesetzt'); xlabel('k'); +figure; stem(kk, abs(X2Hm)); axis([0 N 0 1.1*N/2]); ylabel('|X2_Hm(n)|'); +title('Betragsspektrum bei Hamming-Fensterung'); xlabel('n'); +% ##### EOF ##### diff --git a/common/l3.m b/common/l3.m new file mode 100755 index 0000000..1044ae5 --- /dev/null +++ b/common/l3.m @@ -0,0 +1,29 @@ +% ################################################################################## +% ## Loesung: Reellwertige Exponentialfolge ## +% ## -------------------------------------------------------------------------- ## +% ## Benoetigte(s) m-File(s): lgenexp.m ## +% ################################################################################## + +% ##### Teilaufgabe a: Exponentialfolge erzeugen ##### +[y,k] = lgenexp(0.9, 0, 21); + +figure; stem(k,y); grid; xlabel('k'); ylabel('x(k)'); + +% ##### Teilaufgabe b: Summenformel ueberpruefen ##### +a = 0.9; +k0 = 0; +L = 30; +y = lgenexp(a, k0, L); + +y1 = sum(y); +y2 = (1-a^L)/(1-a); + +% ##### Teilaufgabe c: filter()-Funktion anwenden ##### +a = [1 -0.9]; +b = 1; +d = [1; zeros(20,1)]; +y = filter(b,a,d); +k=0:20; + +figure; stem(k,y); grid; xlabel('k'); ylabel('Amplitude'); +% ##### EOF ##### diff --git a/common/l30.m b/common/l30.m new file mode 100755 index 0000000..1a0d404 --- /dev/null +++ b/common/l30.m @@ -0,0 +1,53 @@ +% ################################################################################## +% ## Loesung: Effizienz des Raderverfahrens ## +% ## -------------------------------------------------------------------------- ## +% ## Benoetigte(s) m-File(s): lrader.m, lsum.m, lakffft.m ## +% ################################################################################## +% Anhand der verwendeten Rechenoperationen (flops) sowie der Rechen- +% zeit wird ein Vergleich der Rechengeschwindigkeit einer AKF- +% Schaetzung nach Rader mit einer direkten Berechnung (Summenformel) +% durchgefuehrt. Die Ergebnisse werden mit diesem m-File aufgetragen. + +N = 2^10; +x = randn(1,N)+j*randn(1,N); +loop = 1:6; + +for l = loop + flops(0); + zeit = cputime; + rxx1 = lrader(2^l,x); + time_rader(l) = cputime-zeit; + flops_rader(l) = flops; + + flops(0); + zeit = cputime; + rxx2 = lsum(2^l,x); + time_sum(l) = cputime-zeit; + flops_sum(l) = flops; + + flops(0); + zeit = cputime; + rxx3 = lakffft(2^l,x); + time_fft(l) = cputime-zeit; + flops_fft(l) = flops; +end; + +disp ('Mean square error (MSE) der AKF-Sch. bzgl. Summenformel:'); +disp (sprintf(' mse_rader: %g', mean(abs(rxx1-rxx2).^2))); +disp (sprintf(' mse_fft : %g', mean(abs(rxx3-rxx2).^2))); + +MM = (2*ones(1,length(loop))).^loop; time_max = max(time_rader); + kflops_max = 300; +figure; plot(MM, time_sum); hold on; plot(MM, time_rader,'--'); +plot(MM, time_fft,':'); axis([0 2^max(loop) 0 1.1*time_max]); xlabel('M'); +title('Zur AKF-Schaetzung benoetigte Rechenzeit'); ylabel('CPU-Zeit in Sek.'); +pos=0.9*time_max; text(38,pos,'___'); text(45,pos,'Summe'); +pos=0.64*time_max; text(38,pos,'- -'); text(45,pos,'Rader'); +pos=0.46*time_max; text(38,pos,'. . .'); text(45,pos,'FFT'); +figure; plot(MM, flops_sum/1000); hold on; plot(MM, flops_rader/1000,'--'); +plot(MM, flops_fft/1000,':'); axis([0 2^max(loop) 0 1.1*300]); +title('Erforderlicher Rechenaufwand in kflops'); ylabel('kflops'); xlabel('M'); +pos=0.22*kflops_max; text(38,pos, '- -'); text(45,pos,'Rader'); +pos=0.95*kflops_max; text(38,pos, '___'); text(45,pos,'Summe'); +pos=0.67*kflops_max; text(38,pos,'. . . '); text(45,pos,'FFT'); +% ##### EOF ##### diff --git a/common/l31.m b/common/l31.m new file mode 100755 index 0000000..f4ee171 --- /dev/null +++ b/common/l31.m @@ -0,0 +1,45 @@ +% ################################################################################## +% ## Loesung: Schaetzung des Leistungsdichtespektrums mit Periodogramm ## +% ## -------------------------------------------------------------------------- ## +% ## Benoetigte(s) m-File(s): lper.m ## +% ################################################################################## +% Teilaufg. a: Lineare Filterung +N = 2^12; % mittelwertfreies, gaussvert., weisses Rauschen +N_vor=500; %Anzahl der Abtastwerte um das System in den eingeschw. + %Zustand zu setzen +n = randn(1,N+N_vor);% der Leistung eins +h1 = [1 -4 6 -4 +1]; x1_all = filter(h1,1,n); x1=x1_all(N_vor+1:length(x1_all)); +h2 = [1 0.8]; x2_all = filter(1,h2,n); x2=x2_all(N_vor+1:length(x2_all)); + +% Teilaufg. b: "Wahres" Leistungsdichtespektrum fuer die beiden Filter +NFFT = 2^12; +Sx1x1 = abs(fft(h1,NFFT)).^2; axis_Sx1=[0 1 0 2*max(Sx1x1)]; +Sx2x2 = ones(1,NFFT)./abs(fft(h2,NFFT)).^2; axis_Sx2=[0 1 0 1.5*max(Sx2x2)]; +fT = 0:1/NFFT:1-1/NFFT; + +figure; plot(fT,Sx1x1); axis([0 1 0 1.1*max(Sx1x1)]);xlabel('Omega/2pi = f T'); +title('Wahres Leistungsdichtespektrum Sx1x1 von x1');ylabel('Sxx(exp(j*Omega))'); +figure; plot(fT,Sx2x2); axis([0 1 0 1.1*max(Sx2x2)]);xlabel('Omega/2pi = f T'); +title('Wahres Leistungsdichtespektrum Sx2x2 von x2');ylabel('Sxx(exp(j*Omega))'); +input('....Fuer naechste Teilaufgabe: RETURN druecken'); + +% Teilaufg. c: Periodogramme und wahre Leistungsdichtespektren zeichnen +fT = 0:1/NFFT:1-1/NFFT; +for k=[2^6, 2^8, 2^10, 2^12] + Sxx=lper(x1(1:k) ,NFFT); + + figure; plot(fT,Sxx); hold on; ylabel('Per(Om)'); + plot(fT,Sx1x1,'--'); xlabel('Omega/2pi = f T'); + title(sprintf('Periodogramm von Sx1x1, N=%d', k)); axis(axis_Sx1); +end; + +input('....Zum Fortfahren: RETURN druecken'); + +for k=[2^6, 2^8, 2^10, 2^12] + Sxx=lper(x2(1:k) ,NFFT); + + figure; plot(fT,Sxx); hold on; ylabel('Per(Om)'); + plot(fT,Sx2x2,'--'); xlabel('Omega/2pi = f T'); + title(sprintf('Periodogramm von Sx2x2, N=%d', k)); axis(axis_Sx2); +end; +% ##### EOF ##### diff --git a/common/l32.m b/common/l32.m new file mode 100755 index 0000000..4ee8aa1 --- /dev/null +++ b/common/l32.m @@ -0,0 +1,53 @@ +% ################################################################################## +% ## Loesung: Schaetzung des Leistungsdichtespektrums nach Blackman-Tuckey ## +% ## -------------------------------------------------------------------------- ## +% ## Benoetigte(s) m-File(s): lblack.m und darin: lrader.m ## +% ################################################################################## + +N = 2^12; +N_vor=500; %Anzahl der Abtastwerte um das System + %in den eingeschw. Zustand zu setzen +n = randn(1,N); % Modellrauschprozesse +h1 = [1 -4 6 -4 +1]; x1_all = filter(h1,1,n); x1=x1_all(N_vor+1:length(x1_all)); +h2 = [1 0.8]; x2_all = filter(1,h2,n); x2=x2_all(N_vor+1:length(x2_all)); + +% "Wahres" Leistungsdichtespektrum zum Vergleich berechnen +NFFT = 2^12; +Sx1x1 = abs(fft(h1,NFFT)).^2; axis_Sx1=[0 1 0 1.35*max(Sx1x1)]; +Sx2x2 = ones(1,NFFT)./abs(fft(h2,NFFT)).^2; axis_Sx2=[0 1 0 1.35*max(Sx2x2)]; +fT = 0:1/NFFT:1-1/NFFT; + +% Schaetzung des Leistungsdichtespektrums fuer M=8, 32, 128 +win = 'triang'; + +Sxx=lblack(x1, 8,win,NFFT); +figure; plot(fT,Sx1x1,'--'); xlabel('Omega/2pi = f T'); +hold on; plot(fT,Sxx); ylabel('Sxx(exp(j*Omega))'); +title('Bl.-Tukey-Sch. von Sx1x1(exp(j*Omega)), M=8 (Dreieck)');axis(axis_Sx1); + +Sxx=lblack(x1, 32,win,NFFT); +figure; plot(fT,Sx1x1,'--'); xlabel('Omega/2pi = f T'); +hold on; plot(fT,Sxx); ylabel('Sxx(exp(j*Omega))'); +title('Bl.-Tukey-Sch. von Sx1x1(exp(j*Omega)), M=32 (Dreieck)'); axis(axis_Sx1); + +Sxx=lblack(x1, 128,win,NFFT); +figure; plot(fT,Sx1x1,'--'); xlabel('Omega/2pi = f T'); +hold on; plot(fT,Sxx); ylabel('Sxx(exp(j*Omega))'); +title('Bl.-Tukey-Sch. von Sx1x1(exp(j*Omega)), M=128 (Dreieck)'); axis(axis_Sx1); +input('....Zum Fortfahren: RETURN druecken'); + +Sxx=lblack(x2, 8,win,NFFT); +figure; plot(fT,Sx2x2,'--'); xlabel('Omega/2pi = f T'); +hold on; plot(fT,Sxx); ylabel('Sxx(exp(j*Omega))'); +title('Bl.-Tukey-Sch. von Sx2x2(exp(j*Omega)), M=8 (Dreieck)');axis(axis_Sx2); + +Sxx=lblack(x2, 32,win,NFFT); +figure; plot(fT,Sx2x2,'--'); xlabel('Omega/2pi = f T'); +hold on; plot(fT,Sxx); ylabel('Sxx(exp(j*Omega))'); +title('Bl.-Tukey-Sch. von Sx2x2(exp(j*Omega)), M=32 (Dreieck)'); axis(axis_Sx2); + +Sxx=lblack(x2, 128,win,NFFT); +figure; plot(fT,Sx2x2,'--'); xlabel('Omega/2pi = f T'); +hold on; plot(fT,Sxx); ylabel('Sxx(exp(j*Omega))'); +title('Bl.-Tukey-Sch. von Sx1x1(exp(j*Omega)), M=128 (Dreieck)'); axis(axis_Sx2); +% ##### EOF ##### diff --git a/common/l33.m b/common/l33.m new file mode 100755 index 0000000..0e1d97a --- /dev/null +++ b/common/l33.m @@ -0,0 +1,52 @@ +% ################################################################################## +% ## Loesung: Schaetzung des Leistungsdichtespektrums nach Bartlett und Welch ## +% ## -------------------------------------------------------------------------- ## +% ## Benoetigte(s) m-File(s): lwelch.m ## +% ################################################################################## + +N = 2^10; N_vor=500; n = randn(1,N+N_vor); % Modellrauschprozesse +h1 = [1 -4 6 -4 +1]; x1_all = filter(h1,1,n); x1=x1_all(N_vor+1:length(x1_all)); +h2 = [1 0.8]; x2_all = filter(1,h2,n); x2=x2_all(N_vor+1:length(x2_all)); +NFFT = 2^10; % "Wahres" Leistungsdichtespektrum zum Vergleich +Sx1x1 = abs(fft(h1,NFFT)).^2; axis_Sx1=[0 1 0 1.5*max(Sx1x1)]; +Sx2x2 = ones(1,NFFT)./abs(fft(h2,NFFT)).^2; axis_Sx2=[0 1 0 1.5*max(Sx2x2)]; +fT = 0:1/NFFT:1-1/NFFT; +v_K = [1 4 16 64]; + +win='boxcar'; % BARTLETT fuer x1 +for k=v_K + S_xx=lwelch(x1,k,win,NFFT); + figure; plot(fT,S_xx); hold on; ylabel('Sxx(exp(j*Omega))'); + plot(fT,Sx1x1,'--'); axis(axis_Sx1); xlabel('Omega/2pi = f T'); + title(sprintf('Bartlett-Sch. von Sx1x1(exp(j*Om.)); K=%d und L=%d/K',k, N)); + text(0.73,900,'== Periodogramm'); +end; +input('......zum Fortfahren: RETURN druecken'); + +win='hamming'; % WELCH-Verf. mit Hamming fuer x1 +for k=v_K + S_xx=lwelch(x1,k,win,NFFT); + figure; plot(fT,S_xx); hold on; ylabel('Sxx(exp(j*Omega))'); + plot(fT,Sx1x1,'--'); axis(axis_Sx1); xlabel('Omega/2pi = f T'); + title(sprintf('Bartlett-Sch. von Sx1x1(exp(j*Om.)); K=%d und L=%d/K',k, N)); +end; +input('......zum Fortfahren: RETURN druecken'); + +win='boxcar'; % BARTLETT fuer x2 +for k=v_K + S_xx=lwelch(x2,k,win,NFFT); + figure; plot(fT,S_xx); hold on; ylabel('Sxx(exp(j*Omega))'); + plot(fT,Sx2x2,'--'); axis(axis_Sx2); xlabel('Omega/2pi = f T'); + title(sprintf('Bartlett-Sch. von Sx2x2(exp(j*Om.)); K=%d und L=%d/K',k, N)); + text(0.73,900,'== Periodogramm'); +end; +input('......zum Fortfahren: RETURN druecken'); + +win='hamming'; % WELCH-Verf. mit Hamming fuer x2 +for k=v_K + S_xx=lwelch(x2,k,win,NFFT); + figure; plot(fT,S_xx); hold on; ylabel('Sxx(exp(j*Omega))'); + plot(fT,Sx2x2,'--'); axis(axis_Sx2); xlabel('Omega/2pi = f T'); + title(sprintf('Bartlett-Sch. von Sx2x2(exp(j*Om.)); K=%d und L=%d/K',k, N)); +end; +% ##### EOF diff --git a/common/l34.m b/common/l34.m new file mode 100755 index 0000000..185dae7 --- /dev/null +++ b/common/l34.m @@ -0,0 +1,90 @@ +% ################################################################################## +% ## Loesung: Yule-Walker und Burg-Algorithmus ## +% ## -------------------------------------------------------------------------- ## +% ## Benoetigte(s) m-File(s): lywex.m, lburg.m ## +% ################################################################################## + +NFFT = 2^10; N = 2^12; +p=[1 8 32]; MA = [1 -4 6 -4 1]; AR = 1; + +% Teilaufg. a +q = length(MA)-1; +disp(sprintf('Wahre AKF (rxx(-%d),...,rxx(0),...,rxx(%d)):', q, q)); +rxx = xcorr(MA).' + +% Teilaufg. b und c: +%%%%%%%%%%%%%%%%%%%%%%%%%% Yule-Walker %%%%%%%%%%%%%%%%%%%%%%%%% +% Ausgabe des tatsaechlichen ARMA- und des "geschaetzten" AR-Betragspektrums +Omega_norm = 0:1/NFFT:1-1/NFFT; +for k=1:3 + [Sxx_ar, Sxx_arma, ar]=lywex(MA,AR,p(k),NFFT); + + figure; plot(Omega_norm, abs(Sxx_arma),'--'); hold on; + plot(Omega_norm, abs(Sxx_ar)); hold off; + axis([0 1 0 1.25*max(Sxx_arma)]); ylabel('Sxx(Omega)'); + title(sprintf('LDS von YW-AR(%d)', p(k))); xlabel('Omega/2pi'); + figure; zplane(MA,ar); title(sprintf('Pole von YW-AR(%d)',p(k))); + xlabel('Realteil'); ylabel('Imaginaerteil'); +end; + +input('....naechste Teilaufgabe: RETURN druecken'); + +% Teilaufg. d: +%%%%%%%%%%%%%%%%%%%%%%%%%% Burg-Methode %%%%%%%%%%%%%%%%%%%%%%%%% +% Ausgabe des tatsaechlichen ARMA- und des geschaetzten AR-Betragspektrums +Omega_norm = 0:1/NFFT:1-1/NFFT; +Sxx_arma = abs(fft(MA,NFFT)./fft(AR,NFFT)).^2; + +for k=1:3 + ar=lburg(N,p(k),MA,AR,NFFT); + + Sxx_ar = abs(ones(1,NFFT)./fft(ar,NFFT)).^2; + Sxx_ar = Sxx_ar/sum(Sxx_ar).*sum(Sxx_arma); + figure; plot(Omega_norm, Sxx_arma,'--'); hold on; + plot(Omega_norm, Sxx_ar); hold off; + axis([0 1 0 1.25*max(Sxx_arma)]); ylabel('Sxx(Omega)'); + title(sprintf('LDS von Burg-AR(%d)', p(k))); xlabel('Omega/2pi'); + figure; zplane(MA,ar); title(sprintf('Pole von Burg-AR(%d)',p(k))); + xlabel('Realteil'); ylabel('Imaginaerteil'); +end; + +input('....naechste Teilaufgabe: RETURN druecken'); + +% Teilaufg. e: +MA=1; AR = [1 0.8]; + +%%%%%%%%%%%%%%%%%%%%%%%%%% Yule-Walker %%%%%%%%%%%%%%%%%%%%%%%%% +% Ausgabe des tatsaechlichen ARMA- und des "geschaetzten" AR-Betragspektrums +Omega_norm = 0:1/NFFT:1-1/NFFT; +for k=1:3 + [Sxx_ar, Sxx_arma, ar]=lywex(MA,AR,p(k),NFFT); + + figure; plot(Omega_norm, abs(Sxx_arma),'--'); hold on; + plot(Omega_norm, abs(Sxx_ar)); hold off; + axis([0 1 0 1.25*max(Sxx_arma)]); ylabel('Sxx(Omega)'); + title(sprintf('LDS von YW-AR(%d)', p(k))); xlabel('Omega/2pi'); + figure; zplane(MA,ar); title(sprintf('Pole von YW-AR(%d)',p(k))); + xlabel('Realteil'); ylabel('Imaginaerteil'); +end; + + +input('....Zum Fortfahren: RETURN druecken'); + +%%%%%%%%%%%%%%%%%%%%%%%%%% Burg-Methode %%%%%%%%%%%%%%%%%%%%%%%%% +% Ausgabe des tatsaechlichen ARMA- und des geschaetzten AR-Betragspektrums +Omega_norm = 0:1/NFFT:1-1/NFFT; +Sxx_arma = abs(fft(MA,NFFT)./fft(AR,NFFT)).^2; + +for k=1:3 + + ar=lburg(N,p(k),MA,AR,NFFT); + Sxx_ar = abs(ones(1,NFFT)./fft(ar,NFFT)).^2; + Sxx_ar = Sxx_ar/sum(Sxx_ar).*sum(Sxx_arma); + figure; plot(Omega_norm, Sxx_arma,'--'); hold on; + plot(Omega_norm, Sxx_ar); hold off; + axis([0 1 0 1.25*max(Sxx_arma)]); ylabel('Sxx(Omega)'); + title(sprintf('LDS von Burg-AR(%d)', p(k))); xlabel('Omega/2pi'); + figure; zplane(AR,ar); title(sprintf('Pole von Burg-AR(%d)',p(k))); + xlabel('Realteil'); ylabel('Imaginaerteil'); +end; +% ##### EOF ##### diff --git a/common/l4.m b/common/l4.m new file mode 100755 index 0000000..d8aa4c3 --- /dev/null +++ b/common/l4.m @@ -0,0 +1,21 @@ +% ################################################################################## +% ## Loesung: Komplexwertige Exponentialfolge ## +% ## -------------------------------------------------------------------------- ## +% ## Benoetigte(s) m-File(s): lgenexp.m ## +% ################################################################################## + +% ##### Teilaufgabe a: komplexe Exponentialfolge erzeugen ##### +z0 = 0.9*( cos(pi/4) + j*sin(pi/4) ); +[x,k] = lgenexp(z0, 0, 21); + +figure; stem(k,real(x)); xlabel('k'); ylabel('Amplitude'); title('Realteil'); grid; +figure; stem(k,imag(x)); xlabel('k'); ylabel('Amplitude'); title('Imaginaerteil'); +grid; + +% ##### Teilaufgabe b: Imaginaerteil ueber Realteil ##### +z0 = 0.9*( cos(pi/4) + j*sin(pi/4) ); +x = lgenexp(z0, 0, 21 ); + +figure; plot(x); grid; xlabel('Realteil'); ylabel('Imaginaerteil'); +title('Phi = 45 [Grad]'); +% ##### EOF ##### diff --git a/common/l5.m b/common/l5.m new file mode 100755 index 0000000..4d3d5c5 --- /dev/null +++ b/common/l5.m @@ -0,0 +1,12 @@ +% ################################################################################## +% ## Loesung: Differenzengleichungen ## +% ################################################################################## + +% ##### Teilaufgabe c: 128 Werte der Impulsantwort ##### +b = [0.3 0.6 0.3]; +a = [1 0 0.9]; +k=[0:127]'; + +impz(b,a,k); %Berechnung der Impulsantwort +title('Impulsantwort'); grid; xlabel('k'); ylabel('Amplitude'); +% ##### EOF ##### diff --git a/common/l6.m b/common/l6.m new file mode 100755 index 0000000..d54743c --- /dev/null +++ b/common/l6.m @@ -0,0 +1,83 @@ +% ################################################################################## +% ## Loesung: Impuls-, Sprungantworten, Frequenzgang ## +% ################################################################################## + +% ##### Teilaufgabe a: Bestimmung der Impulsantwort ##### +x = zeros(111,1); +x(11) = 1; +b = [1 0.5]; +a = [1 -(1.8*cos(pi/16)) 0.81]; +y = filter(b,a,x); +k=[-10:100]'; + +figure; stem(k,y); axis([-10 100 -1 5]); grid; xlabel('k'); +ylabel('Amplitude'); title('Impulsantwort'); + +% ##### Teilaufgabe b: Bestimmung der Sprungantwort ##### +x = [zeros(10,1) ; ones(101,1)]; +b = [1 0.5]; +a = [1 -(1.8*cos(pi/16)) 0.81]; +y = filter(b,a,x); +k=[-10:100]'; + +figure; stem(k,y); axis([-10 100 0 40]); grid; xlabel('k'); +ylabel('Amplitude'); title('Sprungantwort'); + +% ##### Teilaufgabe c: Erste Differenz der Sprungfunktion ##### +d = zeros(111,1); +d(11) = 1; +b = [1 0.5]; +a = [1 -(1.8*cos(pi/16)) 0.81]; +h = filter(b,a,d); +s = [zeros(10,1) ; ones(101,1)]; +g = filter(b,a,s); +h2 = diff(g,1); +k = [-10:100]'; + +figure; stem(k(2:length(k)),h2); grid; xlabel('k'); ylabel('Amplitude'); +title('Erste Differenz der Sprungantwort'); axis([-10 100 -1 5]); +disp('Zum Fortfahren: Druecken Sie eine beliebige Taste...'); pause; + +% ##### Teilaufgabe d: Ausgangssignal ##### +x = exp(j*[1:100]*(pi/3)); +b = [1 0.5]; +a = [1 -(1.8*cos(pi/16)) 0.81]; +y = filter(b,a,x); + +figure; plot(real(y)); grid; xlabel('k'); title('Realteil'); +figure; plot(imag(y)); grid; xlabel('k'); title('Imaginaerteil'); + +% ##### Teilaufgabe e: Uebertragungsfunktionen von LTI Systemen ##### +% ### 2: +b = [1 0.5]; +a = [1 -(1.8*cos(pi/16)) 0.81]; +[H w] = freqz(b,a,128); + +figure; semilogy(w/pi,abs(H)); grid; xlabel('Omega/pi'); +ylabel('| H(exp(j*Omega)) |'); +figure; plot(w/pi,angle(H)); grid; xlabel('Omega/pi'); +ylabel('b(Omega)'); + +% ### 3: +b = [1 0.5]; +a = [1 -(1.8*cos(pi/16)) 0.81]; +[H w] = freqz(b,a,128,'whole'); + +figure; semilogy(w/pi,abs(H)); grid; xlabel('Omega/pi'); +ylabel('| H(exp(j*Omega)) |'); +figure; plot(w/pi,angle(H)); grid; xlabel('Omega/pi'); +ylabel('b(Omega)'); + +% ##### Teilaufgabe f: Pol-Nullstellen-Berechnung ##### +b = [1 0.5]; +a = [1 -(1.8*cos(pi/16)) 0.81]; +Nullstellen = roots(b); +Polstellen = roots(a); + +% ##### Teilaufgabe h: Pol-Nullstellen-Diagramm ##### +b = [1 0.5]; +a = [1 -(1.8*cos(pi/16)) 0.81]; + +figure; zplane(b,a);title('Pol-Nullstellen-Diagramm'); +xlabel('Realteil'); ylabel('Imaginaerteil'); +% ##### EOF ##### diff --git a/common/l7.m b/common/l7.m new file mode 100755 index 0000000..c50aec6 --- /dev/null +++ b/common/l7.m @@ -0,0 +1,65 @@ +% ################################################################################## +% ## Loesung: Fourier-Transformation zeitdiskreter Signale ## +% ## -------------------------------------------------------------------------- ## +% ## Benoetigte(s) m-File(s): ldtft.m ## +% ################################################################################## + +% ##### Teilaufgabe a: Betrag und Phase ##### +a = 0.88 * exp( j * 2*pi/5 ); +k = 0:40; +s = a.^k; +[S W] = ldtft(s, 128); + +figure; plot(W/pi, abs(S) ); grid; title('Betrag'); xlabel('Omega/pi'); +figure; plot(W/pi, angle(S) ); grid; title('Phase'); xlabel('Omega/pi'); +ylabel('Radian'); + +% ##### Teilaufgabe b: Rechteckimpulsfolge ##### +% ### 1: +r =[ones(12,1); zeros(12,1)]'; +k=0:23; + +figure; stem(k,r); grid; xlabel('k'); ylabel('Amplitude'); +title('Rechteckimpulsfolge'); + +% ### 3: +r =[ones(12,1); zeros(12,1)]'; +[R,W]=ldtft(r,256); + +figure; plot(W/pi, abs(R) ); grid; title('Betrag'); xlabel('Omega/pi'); +figure; plot(W/pi, angle(R) ); grid; title('Phase'); xlabel('Omega/pi'); +ylabel('Radian'); + +% ##### Teilaufgabe c: Dreieckimpulsfolge ##### +r =[ones(12,1); zeros(12,1)]'; +d = conv(r,r); +[D,W]=ldtft(d,256); + +figure; plot(W/pi, abs(D) ); grid; title('Betrag'); xlabel('Omega/pi'); +figure; plot(W/pi, angle(D) ); grid; title('Phase'); xlabel('Omega/pi'); +ylabel('Radian'); + +% ##### Teilaufgabe d: Komplexe Modulation ##### +r =[ones(12,1); zeros(12,1)]'; +k = 0:23; +W0 = 2*pi/4; +M = exp( sqrt(-1)*W0*k ); +r = r.*M; +[R,W]=ldtft(r,256); + +figure; plot(W/pi, abs(R) ); grid; title('Betrag'); xlabel('Omega/pi'); +figure; plot(W/pi, angle(R) ); grid; title('Phase'); xlabel('Omega/pi'); +ylabel('Radian'); + +% ##### Teilaufgabe e: Reelle Modulation ##### +r =[ones(12,1); zeros(12,1)]'; +k = 0:23; +W0 = 2*pi/4; +M = exp( sqrt(-1)*W0*k ); +r = r.* real(M); +[R,W]=ldtft(r,256); + +figure; plot(W/pi, abs(R) ); grid; title('Betrag'); xlabel('Omega/pi'); +figure; plot(W/pi, angle(R) ); grid; title('Phase'); xlabel('Omega/pi'); +ylabel('Radian'); +% ##### EOF ##### diff --git a/common/l8.m b/common/l8.m new file mode 100755 index 0000000..1493d19 --- /dev/null +++ b/common/l8.m @@ -0,0 +1,11 @@ +% ################################################################################## +% ## Loesung: Abtastung ## +% ## -------------------------------------------------------------------------- ## +% ## Benoetigte(s) m-File(s): lcsin.m ## +% ################################################################################## + +[y,t]= lcsin(50, 1200, pi/4, 0, 7e-3, 8000); %Zeitdiskretes Sinussignal + +figure; stem(t,y);grid;xlabel('Zeit [sek]');ylabel('Amplitude');title('s_k(kT_A)'); +figure; stem(y); grid;xlabel('k'); ylabel('Amplitude');title('s(k)'); +% ##### EOF ##### diff --git a/common/l9.m b/common/l9.m new file mode 100755 index 0000000..3527330 --- /dev/null +++ b/common/l9.m @@ -0,0 +1,84 @@ +% ################################################################################## +% ## Loesung: Abtastung und Rekonstruktion ## +% ## -------------------------------------------------------------------------- ## +% ## Benoetigte(s) m-File(s): lzerofill.m, lafplot.m, ldtft.m ## +% ################################################################################## + +% ##### Teilaufgabe a: Analoges Signal ##### +fsim = 80000; +delta_t = (0:1/fsim:0.02)'; +f0 = 2000; +s = cos(2*pi*f0*delta_t); + +figure; plot(delta_t,s); xlabel('Zeit (sek)'); ylabel('Amplitude'); +title('Analoges Signal'); + +% ##### Teilaufgabe b: Spektrum des analogen Signals ##### +fsim = 80000; +delta_t = (0:1/fsim:0.02)'; +f0 = 2000; +s = cos(2*pi*f0*delta_t); + +figure; lafplot(s,1/80000); + +% ##### Teilaufgabe c: Abtastung ##### +fsim = 80000; +delta_t = (0:1/fsim:0.02)'; +f0 = 2000; +s = cos(2*pi*f0*delta_t); +fa = 8000; +l = fsim/fa; +index = (1:l:length(s))-1; +index(1) = []; +sd = s(index); + +figure; stem(sd(1:50)); grid; xlabel('k'); ylabel('Amplitude'); +title('Diskretes Signal'); + +% ##### Teilaufgabe d: Spektrum des zeitdisk. Signals ##### +fsim = 80000; +delta_t = (0:1/fsim:0.02)'; +f0 = 2000; +s = cos(2*pi*f0*delta_t); +fa = 8000; +l = fsim/fa; +index = (1:l:length(s))-1; +index(1) = []; +sd = s(index); +[H,W] = ldtft(sd,512); + +figure; plot(W/pi,abs(H)); grid; xlabel('Omega/pi'); ylabel('Amplitude'); +title('Zeitdiskrete Fourier-Transformation (Betrag)'); + +% ##### Teilaufgabe e: Rekonstruktions-TP ##### +fa = 8000; +fsim = 80000; +fcut = fa/fsim; +[b,a] =cheby2(9,60,fcut); + +[H,W]=freqz(b,a,128,fsim); +figure; semilogy(W,abs(H)); grid; xlabel('Frequenz in Hz'); ylabel('Betrag'); +title('Rekonstruktions-Tschebyscheff-TP'); +figure; plot(W,angle(H)); grid; xlabel('Frequenz in Hz'); ylabel('Phase'); +title('Rekonstruktions-Tschebyscheff-TP'); + +% ##### Teilaufgabe f: Analoges Signal ##### +fsim = 80000; % analoges Signal +delta_t = (0:1/fsim:0.02)'; +f0 = 2000; +s = cos(2*pi*f0*delta_t); +fa = 8000; % Abtastung +l = fsim/fa; +index = (1:l:length(s))-1; +index(1) = []; +sd = s(index); % diskretes Signal +sa = lzerofill(sd,l); % Rekonstruktion/Einfuegen von Nullen +fsim = 80000; % Tiefpass entwerfen +fcut = fa/fsim; % siehe Loesung +[b,a]=cheby2(9,60,fcut); +sr = filter(b,a,sa); % Rekonstruktion + +figure; plot(sr); title('Rekonstruierter Sinus'); xlabel... +('Zeit in Abtastwerten'); ylabel('Amplitude'); +figure; lafplot(sr,1/80000); ylabel('Amplitude'); +% ##### EOF ##### diff --git a/common/lafplot.m b/common/lafplot.m new file mode 100755 index 0000000..3e18f8c --- /dev/null +++ b/common/lafplot.m @@ -0,0 +1,20 @@ +% ################################################################################## +% ## Funktion lafplot.m ## +% ################################################################################## +% AFPLOT Stellt die Fourier-Transformierte eines "ANALOGEN" Signals dar. +% Aufruf: lafplot( xa, dt ) +% xa : "ANALOGES" Signal +% dt : Abtastinterval fuer die Simulation von xa(t) + +function lafplot( xa, dt ) + +l = length(xa); +Nfft = round( 2 .^ nextpow2(5*l) ); +Xa = fft(xa, Nfft); +range = 0:(Nfft/4); +ff = range/Nfft/dt; + +plot( ff/1000, abs( Xa(1+range) ) ); +title('Zeitkontinuierliche Fourier-Transformation (Betrag)'); +xlabel('Frequenz (kHz)'); grid +% ##### EOF ##### diff --git a/common/lakffft.m b/common/lakffft.m new file mode 100755 index 0000000..7d74986 --- /dev/null +++ b/common/lakffft.m @@ -0,0 +1,24 @@ +% ################################################################################## +% ## Funktion: lakffft AKF-Schaetzung nach IDFT{|DFT{x}|^2} ## +% ################################################################################## +% +% function rxx = lakffft(M,x) +% +% Nicht erwartungstreue AKF-Schaetzung nach IDFT{|DFT{x}|^2}. +% Der Parameter M entspricht dem Maximalwert fuer lambda. +% Der Vektor x mit den Signalwerten darf eine beliebige Dimension +% annehmen und kann als Spalten- oder Zeilenvektor uebergeben werden. +% Die Ausgabe erfolgt als Spaltenvektor der Dimension M sowie nur fuer +% positive lambda-Werte: rxx = [rxx(0), ..., rxx(M-1)]. +% Fuer negative lambda ist eine konj. gerade Ergaenzung vorzunehmen. + +function rxx = lakffft(M,x) + +x = x(:); % Spaltenvektoren erzeugen +N = length(x); +NFFT = pow2(nextpow2(N+M)); +X = fft(x,NFFT); +Sxx = abs(X.^2); +rxx = ifft(Sxx)./N; +rxx = rxx(1:M); +% ##### EOF ##### diff --git a/common/lband.m b/common/lband.m new file mode 100755 index 0000000..13f6b91 --- /dev/null +++ b/common/lband.m @@ -0,0 +1,18 @@ +% ################################################################################## +% ## Funktion: lband.m; Erzeugung eines reellen Bandpass-Signals ## +% ################################################################################## +% +% function b = lband(N,fm_norm,b_norm); +% +% N Datenwerte eines reellen Bandpass-Signals mit der normierten +% Mittenfrequenz fm_norm=fm/fa und der normierten Bandbreite +% b_norm=b/fa werden erzeugt. + +function b = lband(N,fm_norm,b_norm); + +b = zeros(1,N); +for nu = -5:5 + b = b + (1+nu/10)*cos(2*pi*(fm_norm + b_norm*nu/10)*(1:N)); +end +b = b.'/11; +% ##### EOF ##### diff --git a/common/lblack.m b/common/lblack.m new file mode 100755 index 0000000..3c11c8c --- /dev/null +++ b/common/lblack.m @@ -0,0 +1,24 @@ +% ################################################################################## +% ## Funktion: lblack.m; Berechnung des Blackman-Tuckey-Leistungsdichtespektrums ## +% ## -------------------------------------------------------------------------- ## +% ## Benoetigte(s) m-File(s): lrader.m ## +% ################################################################################## +% +% function Sxx = lblack(x,M,windowtype[,NFFT]); +% +% Schaetzung der spektralen Leistungsdichte nach Blackman-Tukey (1958) +% x Datenfolge (Zeilen- oder Spaltenvektor) +% M Anzahl der verwendeten Autokorrelationswerte +% (sollte Zweierpotenz sein): rxx(-(M-1)), ..., rxx(M-1) +% windowtype Fensterfunktion (String, z.B. 'hamming'), +% [NFFT] FFT-Laenge (optional) + +function Sxx = lblack(x,M,windowtype,NFFT) + +rxx = lrader(M,x); % rxx = [rxx(0); ...; rxx(M-1)] +rxx = [flipud(conj(rxx(2:M))); rxx]; % Konj. gerade Ergaenzung + % => Laenge 2M-1 +eval(['fenster= ', windowtype,'(2*M-1);']); % Fensterfunktion der Laenge +if nargin<4, NFFT=nextpow2(length(x)); end; % 2M-1 berechnen +Sxx = abs(fft(rxx.*(fenster),NFFT)); +% ##### EOF ##### diff --git a/common/lburg.m b/common/lburg.m new file mode 100755 index 0000000..12ee7b6 --- /dev/null +++ b/common/lburg.m @@ -0,0 +1,43 @@ +% ################################################################################## +% ## Funktion: lburg.m; Burg-Algorithmus ## +% ################################################################################## +% +% function [ARcoeff,g] = lburg(N,q,MA,AR[,NFFT,[,N_vor]]); +% +% Spektralschaetzung auf der Grundlage des Burg-Algorithmus' +% Input-Argumente: +% N : Laenge des Datenblockes +% q : Grad des Praediktorfehlerfilters +% AR,MA : Koeffizienten des ARMA-Modellfilters +% [NFFT] : FFT-Laenge +% [N_vor]: Anzahl Abtastwerte fuer eingschw. Zustand +% Output-Argumente: +% ARcoeff: Koeffizienten des Burg-AR-Modells (Zeilenvektor) +% g : PARCOR-Koeffizienten + +function [ARcoeff,g] = lburg(N,q,MA,AR,NFFT,N_vor) + +if nargin<5, NFFT=2^10; N_vor=0; end; +if nargin<6, N_vor=0; end; + +x = randn(1,N+N_vor); +y_all= filter(MA,AR,x); +y = y_all(N_vor+1:length(y_all)); + % Initialisierung +e = y; % Vorwaerts-Praediktorfehler +b = y; % Rueckwaerts-Praediktorfehler + + % Burg-Algorithmus +for ii = 1:q + g(ii) = 2*e(2:N+1-ii)*b(1:N-ii)'/(e(2:N+1-ii)*e(2:N+1-ii)'+... +b(1:N-ii)*b(1:N-ii)'); + e2 = e(2:N+1-ii) - g(ii)*b(1:N-ii); + b = b(1:N-ii)-conj(g(ii))*e(2:N+1-ii); + e = e2; +end; + +ARcoeff = [1 -g(1)]; % Levinson-Durbin-Algorithmus +for i = 2:q + ARcoeff = [ARcoeff 0]-fliplr(conj([ARcoeff 0]))*g(i); +end +% ##### EOF ##### diff --git a/common/lcascade.m b/common/lcascade.m new file mode 100755 index 0000000..b67ca93 --- /dev/null +++ b/common/lcascade.m @@ -0,0 +1,15 @@ +% ################################################################################## +% ## Funktion [B,A]=lcascade(b,a) ## +% ################################################################################## +% Kaskadierung eines durch seine Uebertragungsfunktion gegebenes System in eine +% Kaskadenstruktur aus Systemen zweiter Ordnung. Die Zeilen von B und A beinhalten +% die Koeffizienten des Nenners bzw. des Zaehlers der zugehoerigen Systeme +% zweiter Orndung. + +function [B,A]=lcascade(b,a) + +[z,p,k] = tf2zp(b,a); + sos = zp2sos(z,p,k); +B = sos(:,1:3); +A = sos(:,4:6); +% ##### EOF ##### diff --git a/common/lcheby.m b/common/lcheby.m new file mode 100755 index 0000000..edb6cf2 --- /dev/null +++ b/common/lcheby.m @@ -0,0 +1,41 @@ +% ################################################################################## +% ## Funktion: lcheby.m; Berechnung der kausalen Chebychev-Fensterfunktion ## +% ################################################################################## +% +% function [f_dt,F_dt,a_min_dB] = lcheby(N,Omega_s,delta_Omega) +% +% Berechnung der kausalen Chebychev-Fensterfunktion mittels Ruecktrafo +% aus dem Spektralbereich. Diese Routine liefert im Gegensatz zur +% Matlab-Funktion "chebwin" auch fuer gerade N das richtige Ergebnis. +% Inputs: N : Fensterlaenge im Zeitbereich +% Omega_s : Sperrbereichs-Grenzfrequenz +% delta_Omega: Stuetzstellenabstand im Frequenzbereich +% Outputs: f_dt : Kausale Fensterfunktion (reell) +% F_dt : Spektralfunktion des Fensters (komplex) +% a_min_dB : Erreichte Sperrdaempfung in dB. +% Bsp.: N=34; lcheby(N,0.1*pi,2*pi/N/8) + +function [f_dt,F_dt,a_min_dB] = lcheby(N,Omega_s,delta_Omega) + +Omega = 0:delta_Omega:2*pi-delta_Omega; +F_dt = [cosh((N-1)*acosh(cos(Omega/2)/cos(Omega_s/2)))]; +F_dt = F_dt .*exp(-j*(N-1)*Omega/2); % damit f_dt kausal wird +f_dt = real(ifft(F_dt)); % Imaginaerteile (=1e-14) abschneiden + +fn_dt = f_dt/max(f_dt); % reell, max{fn_dt}=1 +Fn_dt = F_dt/max(abs(F_dt)); % kompl, max{|Fn_dt|}=1 = 0dB +aFn_dt_dB= 20*log10(abs(Fn_dt)); +a_min_dB = 20*log10(cosh((N-1)*acosh(1/cos(Omega_s/2)))); + +figure; fT = Omega/(2*pi); f_sT = Omega_s/(2*pi); y_axis = [-a_min_dB-35 10]; +plot(fT, aFn_dt_dB); hold on; ylabel('|F(exp(j*Om.))/F(exp(j*0))| in dB'); +title(sprintf('|Norm. DT-Spektralfunktion|: N=%d, Om_s/2pi=%.2f,dOm/2pi=1/%d',... + N, f_sT, round(2*pi/delta_Omega))); xlabel('Omega/2pi = f T'); +textstring = sprintf('a_{min} = %.1f dB', a_min_dB); +text(0.22, -0.8*(a_min_dB-5), textstring); axis([0 1 y_axis]); +plot([f_sT f_sT],y_axis,':'); +text(f_sT+0.01,-0.1*(a_min_dB-5),'\Omega_{s}/2\pi=f_{s} T'); +figure; bar(0:length(fn_dt)-1, fn_dt); hold on; plot([N-1 N-1],[0 1.1],':'); +title('Normierte Koeffizienten des Dolph-Tschebyscheff-Fensters'); +ylabel('f(k)'); xlabel('k'); axis([-0.5 2*N-0.5 0 1.1]); +% ##### EOF ##### diff --git a/common/lcoefrnd.m b/common/lcoefrnd.m new file mode 100755 index 0000000..0eb8f1e --- /dev/null +++ b/common/lcoefrnd.m @@ -0,0 +1,14 @@ +% ################################################################################## +% ## Funktion [aq,nfa]=lcoefrnd(a,w) ## +% ################################################################################## +% Quantisierung eines Koeffizientenvektors A auf eine gewuenschte Wortlaenge W. + +function [aq,nfa]=lcoefrnd(a,w) + +f = log(max(abs(a)))/log(2); % Normierung von A +n = 2^ceil(f); % durch n (Potenz von 2), so dass, +an = a/n; % 1 >= an >= -1. +aq = lfxquant(an,w); % Quantisierung von a, so dass, + % 1 > aq >= -1; +nfa = n; % Normierungsfaktor +% ##### EOF ##### diff --git a/common/lcsin.m b/common/lcsin.m new file mode 100755 index 0000000..02602ef --- /dev/null +++ b/common/lcsin.m @@ -0,0 +1,21 @@ +% ################################################################################## +% ## Funktion: lcsin.m; Abtastwerte generieren ## +% ################################################################################## +% +% function [y,t]=lcsin(A,f,phi,t1,t2,Fs) +% +% Generiert Abtastwerte eines zeitkontinuierlichen +% Sinussignals. +% A: Amplitude t1: Startzeit in Sekunden +% f: Frequenz in Hertz t2: Stopzeit in Sekunden +% phi: Initialisierungsphase Fs: Abtatstfrequenz in Hertz + +function [y,t]=lcsin(A,f,phi,t1,t2,Fs) + + if nargin ~= 6 % Ueberpruefen der Eingabeparameter + error('Falsche Anzahl an Eingabeparameter!') +end +t = (t1:1/Fs:t2)'; +y = A*sin(2*pi*f*t+phi); +% ##### EOF ##### + diff --git a/common/ldsin.m b/common/ldsin.m new file mode 100755 index 0000000..977b2ef --- /dev/null +++ b/common/ldsin.m @@ -0,0 +1,14 @@ +% ################################################################################## +% ## Funktion [y,k]=ldsin(A,w,phi,k1,k2) ## +% ################################################################################## +% Generiert ein Sinussignal endlicher Laenge mit Amplitude A, Frequenz W, Phase PHI +% und Zeitindex k1<= k <= k2: y = A*sin(W*k+PHI). + +function [y,k]=ldsin(A,w,phi,k1,k2) + +if nargin ~= 5 % Ueberpruefen der Eingabeparameter + error('Falsche Anzahl an Eingabeparametern!') +end +k = k1:k2; +y = A*sin(w*k+phi); +% ##### EOF ##### diff --git a/common/ldtft.m b/common/ldtft.m new file mode 100755 index 0000000..368c457 --- /dev/null +++ b/common/ldtft.m @@ -0,0 +1,24 @@ +% ################################################################################## +% ## Funktion [H, W] = ldtft(h, N) ## +% ################################################################################## +% DTFT Berechnet die DTFT an N aequidistanten Frequenzstuetzstellen. +% Aufruf: [H, W] = ldtft(h, N) +% h : Eingabevektor der Laenge L +% N : Anzahl der Frequenzen im Intervall [-pi,pi); (N >= L) +% H : DTFT Werte (komplexwertig) +% W : Vektor der Stuetzstellen, an denen die DTFT berechnet wurde + +function [H, W] = ldtft(h, N) + +N = fix(N); +L = length(h); +h = h(:); % nur Vektoren +if( N < L ) + error('DTFT: Anzahl Datenwerte muss kleiner sein als Anzahl Stuetzstellen') +end; +W = (2*pi/N) * [ 0:(N-1) ]'; +mid = ceil(N/2) + 1; +W(mid:N) = W(mid:N) - 2*pi; % <--- Verschiebung [pi,2pi) to [-pi,0) +W = fftshift(W); +H = fftshift( fft( h, N ) ); % <--- Verschiebung der negativen Frequenzen +% ##### EOF ##### diff --git a/common/lfft.m b/common/lfft.m new file mode 100755 index 0000000..7714c12 --- /dev/null +++ b/common/lfft.m @@ -0,0 +1,32 @@ +% ################################################################################## +% ## Funktion: lfft.m; Fouriertransformation einer reellen Folge ## +% ################################################################################## +% +% function X1 = lfft(x1) +% +% Transformation einer reellen Folge x1(k) der Laenge 2N durch eine +% N-Punkte FFT +% N.B.: x1 muss ein Zeilenvektor sein. + +function X1 = lfft(x1) + +lx = length(x1); +if rem(lx,2)==1 % Ungerade Anzahl von Datenwerten + disp('WARNING: Eine gerade Anzahl von Datenwerten ist notwendig:') + disp(' das letzte Datum wird nicht verwendet.') + x1(lx) = []; + lx = lx-1; +end +N = lx/2; % Laenge von reellem x1: 2N +kk = 1:N; +y = x1(2*kk-1); % Gl. (6.1.7a) +z = x1(2*kk); % Gl. (6.1.7b) +x = y + j*z; % Gl. (6.1.8) + % Laenge von komplexem x: N +X = fft(x); % Laenge von X: N +X = [X X(1)]; % Laenge von X: N+1 +Y = (X + conj(fliplr(X)))/2; % Gl. (6.1.5a) +Z = (X - conj(fliplr(X)))/2/j; % Gl. (6.1.5b) +X1 = Y + exp(-j*pi*(0:N)/N).*Z; % Gl. (6.1.12) +X1 = [X1 fliplr(conj(X1(2:N)))]; % Laenge von X1: 2N +% ##### EOF ##### diff --git a/common/lfxquant.m b/common/lfxquant.m new file mode 100755 index 0000000..f7607b9 --- /dev/null +++ b/common/lfxquant.m @@ -0,0 +1,22 @@ +% ################################################################################## +% ## Funktion X=lfxquant(s,bit) ## +% ################################################################################## +% Funktion quantisiert das Signal S auf eine Wortlaenge von BIT bits auf das +% Intervall [-1,1). Dabei wird auf den naechst hoeheren Wert innerhalb des +% Intervalls (Saettigung) aufgerundet. + +function X=lfxquant(s,bit) + +if nargin ~= 2; + error('Aufruf: lfxquant( S, BIT ).'); +end; +if bit <= 0 | abs(rem(bit,1)) > eps; + error('Wortlaenge muss ein positiver Integer sein'); +end; +Plus1 = 2^(bit-1); +X = s * Plus1; +X = round(X); +X = min(Plus1 - 1,X); % Saettigung +X = max(-Plus1,X); +X = X / Plus1; +% ##### EOF ##### diff --git a/common/lge.m b/common/lge.m new file mode 100755 index 0000000..320b5b5 --- /dev/null +++ b/common/lge.m @@ -0,0 +1,18 @@ +%------------------------------------------------------------------- +% +% [l]=lge(x); +% +% alias fuer length(x): l=max(size(x)); +% +%------------------------------------------------------------------- + +% +% author: Markus Hofbauer, +% ISI, ETH Zuerich (Switzerland) +% +% created: 7/2000 +% +%------------------------------------------------------------------- + +function [l]=lge(x); +l=max(size(x)); diff --git a/common/lgenexp.m b/common/lgenexp.m new file mode 100755 index 0000000..995a3ad --- /dev/null +++ b/common/lgenexp.m @@ -0,0 +1,26 @@ +% ################################################################################## +% ## Funktion: lgenexp.m; Exponentialfolge ## +% ################################################################################## +% +% function [y,k]=lgenexp(a,k0,l) +% +% Generiert eine Exponentialfolge b^k +% a: Eingabeskalar +% k0: Startindex +% l: Laenge der Exponentialfolge + +function [y,k]=lgenexp(a,k0,l) + +if nargin ~= 3 + error('Aufruf ist: Y = GENEXP( B, K0, L)'); +end; +if l<=0 + error('Signallaenge muss positiv sein'); +end; +if abs(a)>1 + error('|a| muss kleiner Eins sein'); +end; +k = k0 + [1:l]' - 1; +y = a .^ k; +% ##### EOF ##### + diff --git a/common/lms.m b/common/lms.m new file mode 100755 index 0000000..8699e36 --- /dev/null +++ b/common/lms.m @@ -0,0 +1,48 @@ +% LMS-Algorithmus +% --------------- +% function [ys,e,wf,fpo] = lms(x,d,N,mu,wi) +% +% Parameter: +% x : Eingangssignal +% d : erwünschtes Signal (Referenzsignal) +% N : Filterordnung +% mu : konstante Schrittweite +% wi : Startwerte der Filterkoeffizienten im Zeitbereich +% +% Rückgabewerte: +% ys : Schätzung des erwünschten Signals d aus x +% e : Fehler d - y +% wf : Filterkoeffizienten am Ende der Adaption im Zeitbereich +% fpo : Anzahl benötigten der Floating-Point Operationen + +% ------------------------------------------------------------------------- +% Datum : 20.7.2001 +% Autor : Jens Ahrensfeld +% Thema : Diplomarbeit +% Datei : lms.m +% Benötigte Dateien: keine +% +% ------------------------------------------------------------------------- +function [ys,e,wf,fpo] = lms(x,d,N,mu,wi) + +% ------------------------------------------------------------------------- +% Initialisierung LMS +% ------------------------------------------------------------------------- +wf = wi; +y = zeros(length(x),1); +e = zeros(length(x),1); +flops(0); + +% ------------------------------------------------------------------------- +% LMS Adaptation +% ------------------------------------------------------------------------- +for n = N : length(x) + xs = x(n:-1:n-N+1); + ys(n) = xs'*wf; + e(n) = d(n) - ys(n); + wf = wf + mu*e(n)*xs; +end ; +fpo = flops; + +% ------------------------------------------------------------------------- +% Ende lms.m diff --git a/common/lms1.m b/common/lms1.m new file mode 100755 index 0000000..f057548 --- /dev/null +++ b/common/lms1.m @@ -0,0 +1,53 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% LMS Algorithm % +% % +% Written By: Sundar Sankaran and A. A. (Louis) Beex % +% DSP Research Laboratory % +% Dept. of Electrical and Comp. Engg % +% Virginia Tech % +% Blacksburg VA 24061-0111 % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + +randn('seed', 0) ; +rand('seed', 0) ; + +NoOfData = 8000 ; % Set no of data points used for training +Order = 32 ; % Set the adaptive filter order + +Mu = 0.01 ; % Set the step-size constant + +x = randn(NoOfData, 1) ;% Input assumed to be white +h = rand(Order, 1) ; % System picked randomly +d = filter(h, 1, x) ; % Generate output (desired signal) + +% Initialize LMS + +w = zeros(Order,1) ; + +% LMS Adaptation + +for n = Order : NoOfData + D = x(n:-1:n-Order+1) ; + d_hat(n) = w'*D ; + e(n) = d(n) - d_hat(n) ; + w = w + Mu*e(n)*D ; + w_err(n) = norm(h - w) ; +end ; + +% Plot results + +figure ; +plot(20*log10(abs(e))) ; +title('Learning Curve') ; +xlabel('Iteration Number') ; +ylabel('Output Estimation Error in dB') ; + +figure ; +semilogy(w_err) ; +title('Weight Estimation Error') ; +xlabel('Iteration Number') ; +ylabel('Weight Error in dB') ; diff --git a/common/lms2.m b/common/lms2.m new file mode 100755 index 0000000..47bad0b --- /dev/null +++ b/common/lms2.m @@ -0,0 +1,67 @@ +%=================================================================== +%------------------------------------------------------------------- +% +% Adaptive Filter +% +% Simulationen +% +% [E,W,w,DD]=lms(mu, N, X, D, w_start) +% +% LMS-Algorithmus +% +% E = Fehlersignal E[.] +% W = Filterkoeffizienten im zeitlichen Verlauf +% w = Filterkoeffizienten am Ende der Adaption +% DD = Filterausgang, Schaetzung von D +% mu = Schrittweite +% N = Anz. Filterkoeffizienten = Filterordnung+1 +% X = Filtereingang X[.] +% D = erwuenschtes Signal D[.] +% w_start = Startwerte Filterkoeffizienten +% +%------------------------------------------------------------------- + +% +% author: Peter Wellig, Martin Haenggi, Markus Hofbauer +% ISI, ETH Zuerich (Switzerland) +% +% created: 7/2000 +% +% +%------------------------------------------------------------------- +% +% File : lms.m +% +% Startfile: simX.m +% +%------------------------------------------------------------------- +%=================================================================== + +function [E,W,w,DD]=lms2(mu, N, X, D, w_start) + +% Initialisierungen +adaptlen=length(X); +w = w_start; +W=zeros(N,adaptlen); +DD=zeros(adaptlen,1); +E=zeros(adaptlen,1); + + +% Loop +flops(0); +for i=N:adaptlen + W(:,i) = w; + + x = X(i:-1:i-N+1); % Eingangsvektor (x[k],x[k-1],..,x[k-N+1]) + y = x'*w; % Filterausgang + e = D(i)-y; % Fehler + w = w+mu*e*x; % Aufdatierung der Filterkoeffizienten + + E(i) = e; + DD(i) = y; +end; +flops + + + + diff --git a/common/lper.m b/common/lper.m new file mode 100755 index 0000000..55828c3 --- /dev/null +++ b/common/lper.m @@ -0,0 +1,17 @@ +% ################################################################################## +% ## Funktion: lper.m; Berechnung des Periodogramms ## +% ################################################################################## +% +% function Sxx = lper(x[,NFFT]) +% +% Schaetzung der spektralen Leistungsdichte mittels Periodogramm +% x Datenfolge (Zeilen- oder Spaltenvektor) +% [NFFT] FFT-Laenge (optional) + +function Sxx = lper(x,NFFT) + +N = length(x); +if nargin<2, NFFT=pow2(nextpow2(N)); end; +X = fft(x,NFFT); +Sxx = (abs(X).^2)./N; +% ##### EOF ##### diff --git a/common/lrader.m b/common/lrader.m new file mode 100755 index 0000000..399f962 --- /dev/null +++ b/common/lrader.m @@ -0,0 +1,44 @@ +% ################################################################################## +% ## Funktion: lrader.m ; AKF-Schaetzung mit Raderverfahren ## +% ################################################################################## +% +% function rxx = lrader(M,x) +% +% Nicht erwartungstreue AKF-Schaetzung nach dem Rader-Verfahren laut +% Kammeyer/Kroschel "Digitale Signalverarbeitung", S. 240-247. +% Der Parameter M muss eine Zweierpotenz sein und entspricht dem +% (Maximalwert+1) fuer lambda. Der Vektor x mit den Signalwerten +% darf eine beliebige Dimension annehmen und kann als Spalten- oder +% Zeilenvektor uebergeben werden. Die Ausgabe erfolgt als Spalten- +% vektor der Dimension M sowie nur fuer positive lambda-Werte: +% rxx = [rxx(0), ..., rxx(M-1)]. +% Fuer negative lambda ist eine konj. gerade Ergaenzung vorzunehmen. + +function rxx = lrader(M,x) +x = x(:); % Spaltenvektoren erzeugen +N = length(x); +u = zeros(2*M,1); +u(1:M) = x(1:M); +U = fft(u); +K = floor(N/M); +v = zeros(2*M,1); +Sxx = zeros(2*M,1); +pm = (-1).^(0:(2*M-1)).'; + +for ii = 1:(K-1) % Iteration und Akkumulation + v(1:M) = x((ii*M+1):((ii+1)*M)); + V = fft(v); + Sxx = Sxx + (conj(U) .*(U + (pm .*V))); + U = V; +end; +if K < N/M + b = N - K*M; + v = zeros(2*M,1); + v(1:b) = x((K*M+1):N); + V = fft(v); + Sxx = Sxx + (conj(U) .*(U + (pm .*V))); +end; +Sxx = Sxx + abs(V.^2); +rxx = ifft(Sxx)/N; +rxx = rxx(1:M); +% ##### EOF ##### diff --git a/common/lrader2.m b/common/lrader2.m new file mode 100755 index 0000000..bee8b0c --- /dev/null +++ b/common/lrader2.m @@ -0,0 +1,50 @@ +% ################################################################################## +% ## Funktion: lrader.m ; AKF-Schaetzung mit Raderverfahren ## +% ################################################################################## +% +% function rxx = lrader(M,x) +% +% Nicht erwartungstreue AKF-Schaetzung nach dem Rader-Verfahren laut +% Kammeyer/Kroschel "Digitale Signalverarbeitung", S. 240-247. +% Der Parameter M muss eine Zweierpotenz sein und entspricht dem +% (Maximalwert+1) fuer lambda. Der Vektor x mit den Signalwerten +% darf eine beliebige Dimension annehmen und kann als Spalten- oder +% Zeilenvektor uebergeben werden. Die Ausgabe erfolgt als Spalten- +% vektor der Dimension M sowie nur fuer positive lambda-Werte: +% rxx = [rxx(0), ..., rxx(M-1)]. +% Fuer negative lambda ist eine konj. gerade Ergaenzung vorzunehmen. + +function rxx = lrader2(M,x,y) +x = x(:); % Spaltenvektoren erzeugen +y = y(:); % Spaltenvektoren erzeugen +N = length(x); +u = zeros(2*M,1); +u(1:M) = x(1:M); +U = fft(u); +w = zeros(2*M,1); +w(1:M) = y(1:M); +W = fft(w); +K = floor(N/M); +v = zeros(2*M,1); +Sxx = zeros(2*M,1); +pm = (-1).^(0:(2*M-1)).'; + +for ii = 1:(K-1) % Iteration und Akkumulation + v(1:M) = y((ii*M+1):((ii+1)*M)); + V = fft(v); + Sxx = Sxx + (conj(U) .*(W + (pm .*V))); + W = V; + u(1:M) = x((ii*M+1):((ii+1)*M)); + U = fft(u); +end; +if K < N/M + b = N - K*M; + v = zeros(2*M,1); + v(1:b) = y((K*M+1):N); + V = fft(v); + Sxx = Sxx + (conj(U) .*(W + (pm .*V))); +end; +Sxx = Sxx + 0*abs(V.^2); +rxx = ifft(Sxx)/N; +rxx = rxx(1:M); +% ##### EOF ##### diff --git a/common/lsum.m b/common/lsum.m new file mode 100755 index 0000000..a895b40 --- /dev/null +++ b/common/lsum.m @@ -0,0 +1,26 @@ +% ################################################################################## +% ## Funktion: lsum.m; AKF-Schaetzung mit Summenformel ## +% ################################################################################## +% +% function rxx = lsum(M,x) +% +% Nicht erwartungstreue AKF-Schaetzung durch Vektormultiplikation +% Der Parameter M entspricht dem (Maximalwert-1) fuer lambda. +% Der Vektor x mit den Signalwerten darf eine beliebige Dimension +% annehmen und kann als Spalten- oder Zeilenvektor uebergeben werden. +% Die Ausgabe erfolgt als Spaltenvektor der Dimension M sowie nur fuer +% positive lambda-Werte: rxx = [rxx(0), ..., rxx(M-1)]. +% Fuer negative lambda ist eine konj. gerade Ergaenzung vorzunehmen. + +function rxx = lsum(M,x) + +x = x(:); % Spaltenvektoren erzeugen +N = length(x); +r0 = x'*x; +r = zeros(M-1,1); +for ii = 1:M-1 + r(ii) = [zeros(1,ii) x'] * [x; zeros(ii,1)]; +end +rxx=[r0; r]./N; % Nicht erwartungstreue AKF-Schaetzung +% rxx=[r0; r]./(N:-1:N-(M-1)); % Erwartungstreue AKF-Schaetzung +% ##### EOF ##### diff --git a/common/lwelch.m b/common/lwelch.m new file mode 100755 index 0000000..83510bf --- /dev/null +++ b/common/lwelch.m @@ -0,0 +1,45 @@ +% ################################################################################## +% ## Funktion: lwelch.m; Welch- bzw. Bartlettverfahren ## +% ################################################################################## +% +% function Sxx = lwelch(x,K,windowtype[,NFFT[,overlap]]) +% +% Schaetzung der spektralen Leistungsdichte nach dem Welch-Ansatz +% x Datenfolge (Zeilenvektor) +% K Anzahl der zu mittelnden Periodogramme = #Teilfolgen +% windowtype Fensterfunktion (z.B. 'hamming'=> Welch-Verfahren, +% 'boxcar' => Bartlett-Verfahren). +% [NFFT] FFT-Laenge (optional) +% [overlap] zusaetzliche Mittelung mit den Periodogrammen der +% um fix(L/2) verschobenen Teilfolgen (falls das +% optionale Argument vorhanden ist) +% Beispiel: Sxx = lwelch(x,4,'boxcar',1024); + +function Sxx = lwelch(x,K,windowtype,NFFT,overlap) + +L = floor(length(x)/K); +if nargin<4, NFFT=nextpow2(L); end; + % Fensterfunktion der Laenge +eval(['fenster =',windowtype,'(L);']); % L berechnen (Spaltenvektor) +AL = fenster.'*fenster; % Normierungskonstante +Sxx = zeros(1,NFFT); +kk = (1:L); + +for i = 0:K-1 % Addition der Periodogramme + yi = fenster.' .* x(kk+i*L); % der nicht-ueberlappenden + Yi = fft(yi,NFFT); % Teilfolgen von x + Sxx = Sxx + abs(Yi).^2; +end + +if nargin > 4 % zusaetzliche Addition der + kk = (1:L) + fix(L/2); % Periodogramme der ueber- + for i = 0:K-2 % lappenden Teilfolgen von x + yi = fenster' .* x(kk+i*L); + Yi = fft(yi,NFFT); + Sxx = Sxx + abs(Yi).^2; + end + Sxx = Sxx/AL/(2*K-1); % Normierung auf AL und +else % Anzahl der Mittelungen (=2K-1) + Sxx = Sxx/AL/K; % Normierung auf AL und +end % Anzahl der Mittelungen (=K) +% ##### EOF ##### diff --git a/common/lywex.m b/common/lywex.m new file mode 100755 index 0000000..cd5c5e4 --- /dev/null +++ b/common/lywex.m @@ -0,0 +1,42 @@ +% ################################################################################## +% ## Funktion: lywex.m; AR-(p)-Approximation eines ARMA-Modelles nach Yule-Walker +% ################################################################################## +% +% function [ARcoeff] = lywex(MA,AR,p[,NFFT]) +% +% Berechnung der Koeffizienten eines AR-Modelles aus den exakten +% Autokorrelationskoeffizienten +% Input-Argumente: +% MA,AR : Parameter des ARMA-Modells +% p : Grad des zu bestimmenden AR-Modells +% [NFFT] : FFT-Laenge +% Output-Argument: +% ARcoeff : Koeffizienten des Yule-Walker-AR-Modells (Zeilenvektor) +% Sxx_ar : geschaetztes AR-Betragspektrum +% Sxx_arma : wahres LDS des ARMA-Modells + +function [Sxx_ar, Sxx_arma, ARcoeff] = lywex(MA,AR,p,NFFT) + +if nargin<4, NFFT = 2^10; end; +delta = [1 zeros(1,NFFT-1)]; +h_arma = filter(MA,AR,delta); +H_arma = fft(h_arma); +Sxx_arma = abs(H_arma).^2; % LDS des ARMA-Modells + % Wahre AKF des ARMA-Modells: +r0_arma = h_arma*h_arma'; % rxx(0) +for ii = 1:p % rxx(1)..rxx(p) + r_arma(ii) = [zeros(1,ii) h_arma]*[h_arma zeros(1,ii)]'; +end + +R_arma = toeplitz([r0_arma r_arma(1:length(r_arma)-1)]); +ARcoeff = -inv(R_arma)*r_arma'; +ARcoeff = [1 ARcoeff.']; +h_ar = filter(1, ARcoeff, delta); +H_ar = fft(h_ar); +if p == 0 + sigmak2 = 1; +else + sigmak2 = R_arma(1,1)-r_arma*inv(R_arma)*r_arma'; +end +Sxx_ar = sigmak2 * abs(H_ar).^2; %geschaetztes AR-Betragspektrum +% ##### EOF ##### diff --git a/common/lzerofill.m b/common/lzerofill.m new file mode 100755 index 0000000..b353827 --- /dev/null +++ b/common/lzerofill.m @@ -0,0 +1,24 @@ +% ################################################################################## +% ## Funktion: lzerofill.m; Sequenz mit Nullen fuellen ## +% ################################################################################## +% +% function y = lzerofill(x,l) +% +% Die Funktion fuellt eine Sequenz X mit l-1 Nullen zwischen jedem Wert von X auf. +% x : Eingangssignalvektor +% l : l-1 Nullen zwischen jedem Abtastwert +% y : Ausgangssignalvektor ==> Laenge(y) = l*Laenge(x) +% NOTE: Falls X eine Matrix ist, wird zerofill auf +% jede Spalte angewendet. + +function y = lzerofill(x,l) + +[M,N] = size(x); +if M==1 % <--- fuer Vektoren + y = zeros(1, N*l); + y(1,1:l:N*l) = x; +else % <--- Fuer Matrizen + y = zeros(M*l, N); + y(1:l:M*l,:) = x; +end +% ##### EOF ##### diff --git a/common/lzp1.m b/common/lzp1.m new file mode 100755 index 0000000..e22c0fd --- /dev/null +++ b/common/lzp1.m @@ -0,0 +1,10 @@ +% ################################################################################## +% ## Funktion v = zp1(M,N) ## +% ################################################################################## +% Die Funktion erzeugt eine M*N-Matrix mit Zufallswerten + +function v = zp1(M,N) + +Mc = ones(M,1)*cos((1:N)/2); +v = (rand(M,N)-0.5) .* Mc; +% ##### EOF ##### diff --git a/common/mx.m b/common/mx.m new file mode 100755 index 0000000..f1d9646 --- /dev/null +++ b/common/mx.m @@ -0,0 +1,36 @@ +% Magnitude Estimation von x +% Schätzung des mittleren Betrages von x +% mit kombinierter Anstiegs- und Abfallzeit trf +% --------------------------------------------- +% +% Aufruf: +% function [mx,mxf] = mx(x,trf,mxi) +% +% Parameter: +% x : Eingangsvektor +% trf : Zeitkonstante Anstiegs- und Abfallzeit +% mxi : Initial Condition Anfangsbetrag von x +% +% Rückgabewerte: +% mx : Ausgang mittlerer Betrag von x +% mxf : Endbetrag von x + +% ------------------------------------------------------------------------- +% Datum : 26.07.2001 +% Autor : Jens Ahrensfeld +% Thema : Diplomarbeit +% Datei : mx.m +% Benötigte Dateien: filter, lge, abs +% +% ------------------------------------------------------------------------- +function [mx,mxf] = mx(x,trf,mxi) + +N = lge(x); +x(1) = x(1) + mxi*(1/trf-1); +%x = abs(x); +a = [1,-(1-trf)]; +mx = filter(trf,a,x); +mxf = mx(N); + +% ------------------------------------------------------------------------- +% Ende mx.m diff --git a/common/mx2.m b/common/mx2.m new file mode 100755 index 0000000..901d55f --- /dev/null +++ b/common/mx2.m @@ -0,0 +1,42 @@ +% Magnitude Estimation von x +% Schätzung des mittleren Betrages von x +% mit getrennter Anstiegs- und Abfallzeit tr, tf +% ---------------------------------------------- +% +% Aufruf: +% function [mx,mxf] = mx2(x,tr,tf,mxi) +% +% Parameter: +% x : Eingangsvektor +% tr : Zeitkonstante Anstiegszeit +% tf : Zeitkonstante Abfallzeit +% mxi : Initial Condition Anfangsbetrag von x +% +% Rückgabewerte: +% mx : Ausgang mittlerer Betrag von x +% mxf : Endbetrag von x + +% ------------------------------------------------------------------------- +% Datum : 27.03.2002 +% Autor : Jens Ahrensfeld +% Thema : Diplomarbeit +% Datei : mx2.m +% +% ------------------------------------------------------------------------- +function [mx,mxf] = mx2(x,tr,tf,mxi) + +N = lge(x); +x = abs(x); +mx(1) = mxi; + +for n = 2:N, + if(x(n) > mx(n-1)) + mx(n) = tr*x(n) + (1-tr)*mx(n-1); + else + mx(n) = tf*x(n) + (1-tf)*mx(n-1); + end; +end; +mxf = mx(N); + +% ------------------------------------------------------------------------- +% Ende mx2.m diff --git a/common/newfunc.m b/common/newfunc.m new file mode 100755 index 0000000..e881feb --- /dev/null +++ b/common/newfunc.m @@ -0,0 +1,18 @@ +% Neue Function +% ------------------------------------ +% +% Aufruf: +% Function []=NeueFunktion() +% +% Parameter: +% +% Rückgabewerte: + +% ------------------------------------------------------------------------- +% Datum : ..2001 +% Autor : Jens Ahrensfeld +% Thema : Diplomarbeit +% Datei : newfunc.m +% Benötigte Dateien: +% +% ------------------------------------------------------------------------- diff --git a/common/nlms.m b/common/nlms.m new file mode 100755 index 0000000..f8be27c --- /dev/null +++ b/common/nlms.m @@ -0,0 +1,48 @@ +% NLMS-Algorithmus +% --------------- +% function [ys,e,wf,fpo] = nlms(x,d,N,mu,wi) +% +% Parameter: +% x : Eingangssignal +% d : erwünschtes Signal (Referenzsignal) +% N : Filterordnung +% mu : konstante Schrittweite +% wi : Startwerte der Filterkoeffizienten im Zeitbereich +% +% Rückgabewerte: +% ys : Schätzung des Referenzsignals d aus x +% e : Fehler d - ys +% wf : Filterkoeffizienten am Ende der Adaption im Zeitbereich +% fpo : Anzahl der benötigten Floating-Point Operationen + +% ------------------------------------------------------------------------- +% Datum : 20.7.2001 +% Autor : Jens Ahrensfeld +% Thema : Diplomarbeit +% Datei : nlms.m +% Benötigte Dateien: keine +% +% ------------------------------------------------------------------------- +function [ys,e,wf,fpo] = nlms(x,d,N,mu,wi) + +% ------------------------------------------------------------------------- +% Initialisierung NLMS +% ------------------------------------------------------------------------- +wf = wi; +ys = zeros(length(x),1); +e = zeros(length(x),1); +flops(0); + +% ------------------------------------------------------------------------- +% NLMS Adaptation +% ------------------------------------------------------------------------- +for n = N : length(x) + xs = x(n:-1:n-N+1); + ys(n) = xs'*wf; + e(n) = d(n) - ys(n); + wf = wf + mu*e(n)*xs/(xs'*xs+0.001); +end ; +fpo = flops; + +% ------------------------------------------------------------------------- +% Ende nlms.m diff --git a/common/nlms1.m b/common/nlms1.m new file mode 100755 index 0000000..68c707a --- /dev/null +++ b/common/nlms1.m @@ -0,0 +1,51 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% Normalized LMS Algorithm % +% % +% Written By: Sundar Sankaran and A. A. (Louis) Beex % +% DSP Research Laboratory % +% Dept. of Electrical and Comp. Engg % +% Virginia Tech % +% Blacksburg VA 24061-0111 % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +randn('seed', 0) ; +rand('seed', 0) ; + +NoOfData = 8000 ; % Set no of data points used for training +Order = 32 ; % Set the adaptive filter order + +Mu = 1.0 ; % Set the step-size constant + +x = randn(NoOfData, 1) ;% Input assumed to be white +h = rand(Order, 1) ; % System picked randomly +d = filter(h, 1, x) ; % Generate output (desired signal) + +% Initialize NLMS + +w = zeros(Order,1) ; + +% NLMS Adaptation + +for n = Order : NoOfData + D = x(n:-1:n-Order+1) ; + d_hat(n) = w'*D ; + e(n) = d(n) - d_hat(n) ; + w = w + Mu*e(n)*D/(D'*D) ; + w_err(n) = norm(h - w) ; +end ; + +% Plot results + +figure ; +plot(20*log10(abs(e))) ; +title('Learning Curve') ; +xlabel('Iteration Number') ; +ylabel('Output Estimation Error in dB') ; + +figure ; +semilogy(w_err) ; +title('Weight Estimation Error') ; +xlabel('Iteration Number') ; +ylabel('Weight Error in dB') ; diff --git a/common/pfftcnv3.m b/common/pfftcnv3.m new file mode 100755 index 0000000..fcd5a02 --- /dev/null +++ b/common/pfftcnv3.m @@ -0,0 +1,86 @@ +% pfftcnv2(x,b) +function y = pfftcnv2(x,h,P,S,L) + +N = length(h); +Lx = length(x); +C = 2^ceil(log(L+S*L-1)/log(2)); +Np = S*L; +K = Lx/L; + +y = zeros(1,Lx); +hp = zeros(P,C); +Hp = zeros(P,C); +xk = zeros(1,C); +Xk = zeros(1,C); +Xkps = zeros(S*P,C); +yk = zeros(1,C); +Y = zeros(1,C); +xzp = zeros(1,Lx+C-L); +pidTbl = [2:P*S]; +pidTbl(P*S) = 1; +pStbl = [0:P*S-1]; +pStbl(1) = 10; + +pid = 1; + +% h(k) => hp(k) +for p=1:P, + SLp = (p-1)*Np; + hp(p,1:Np) = h(SLp+1:SLp+Np); +end; + +% hp(k) => Hp(k) +for p=1:P, + Hp(p,1:C) = fft(hp(p,1:C),C); +end; + +xzp(C-L+1:Lx+C-L) = x; + +tic; + +for k=1:K, + + kL = (k-1)*L; + xk(1:C) = xzp(kL+1:kL+C); + + Xkps(pid,1:C) = fft(xk,C); + + pS = pid; + Y = Hp(1,1:C).*Xkps(pid,1:C); + for p=2:P, + for i = 1:S + pS = pStbl(pS); + end; + Y = Y + Hp(p,1:C).*Xkps(pS,1:C); + end; + + yk = ifft(Y,C); + + y(kL+1:kL+L) = yk(C-L+1:C); + + pid = pidTbl(pid); + +end; +toc; + +SUBPLOT(2,1,1), + plot(1:N,h); + n = [1:N]; + XT = [Np*n]; + set(gca,'XTick',XT); + grid; + PlotTitle = sprintf('Filterkoeffizienten h(n), N=%d, P=%d, S=%d, NP=%d',N,P,S,Np); + title(PlotTitle); + xlabel('n'); + ylabel('h'); + +SUBPLOT(2,1,2), + plot(1:K*L,real(y(1:K*L))); + n = [1:N]; + XT = [L*n]; + set(gca,'XTick',XT); + grid; + PlotTitle = sprintf('Filterausgang y(n), K=%d, L=%d, C=%d',K,L,C); + title(PlotTitle); + xlabel('n'); + ylabel('y'); diff --git a/common/pfftfilt.m b/common/pfftfilt.m new file mode 100755 index 0000000..f0c5775 --- /dev/null +++ b/common/pfftfilt.m @@ -0,0 +1,51 @@ +% pfftfilt(x,b) +function y = pfftfilt(x,h,P,S,L) + +N = length(h); +Lx = length(x); +C = 2^ceil(log(L+S*L-1)/log(2)); +Np = S*L; +K = Lx/L; + +y = zeros(Lx,1); +Hp = zeros(C,P); +Xkps = zeros(C,S*P); +yk = zeros(C,1); +Y = zeros(C,1); +xzp = zeros(Lx+C-L,1); +pidTbl = [2:P*S]; +pidTbl(P*S) = 1; +pStbl = [0:P*S-1]; +pStbl(1) = 10; + +pid = 1; + +% h(k) => Hp(k) +for p=1:P, + Hp(1:C,p) = fft([h((p-1)*Np+1:p*Np);zeros(C-Np,1)],C); +end; + +xzp(C-L+1:Lx+C-L) = x; + +for k=1:K, + + kL = (k-1)*L; + + Xkps(1:C,pid) = fft(xzp(kL+1:kL+C),C); + + pS = pid; + Y = Hp(1:C,1).*Xkps(1:C,pid); + for p=2:P, + for i = 1:S + pS = pStbl(pS); + end; + Y = Y + Hp(1:C,p).*Xkps(1:C,pS); + end; + + yk = ifft(Y,C); + + y(kL+1:kL+L) = yk(C-L+1:C); + + pid = pidTbl(pid); + +end; diff --git a/common/pflms.m b/common/pflms.m new file mode 100755 index 0000000..be0ab6d --- /dev/null +++ b/common/pflms.m @@ -0,0 +1,210 @@ +% Partitioned Frequency LMS Adaptive Filter -PFLMS- +% ------------------------------------------------- +% Aufruf: +% Function [dd,e,wf,fpo]=pflms(x,d,mu,N,C,P,S,L,wi) +% +% Parameter: +% x : Eingangsvektor +% d : Referenzvektor +% mu : zeitabhängiger Schrittweitenvektor 0 < mu <= 1.0 +% N : Anzahl Filterkoeffizienten (Hinweis: N = P*S*L) +% C : FFT-Länge, idealerweise eine Zweier-Potenz (0=Auto) +% P : Anzahl der Filterpartitionen +% S : Anzahl der Filtersegmente +% L : Länge eines Verarbeitungsblocks (idealerweise C=L+N/P=2N/P) +% wi : Startwerte der Filterkoeffizienten im Zeitbereich +% +% Rückgabewerte: +% e : Adaptitionsfehler +% dd : Filterausgang +% wf : Filterkoeffizienten am Ende der Adaption im Zeitbereich +% fpo : Anzahl benötigten der Floating-Point Operationen + +% ------------------------------------------------------------------------- +% Datum : 27.3.2002 +% Autor : Jens Ahrensfeld +% Thema : Diplomarbeit +% Datei : pflms.m +% +% ------------------------------------------------------------------------- +function [dd,e,wf,fpo]=pflms(x,d,mu,N,C,P,S,L,wi) + +% ------------------------------------------------------------------------- +% Überprüfen der Parameter +% ------------------------------------------------------------------------- +if (N ~= fix(N)), + error('Fehler: N muss eine ganze Zahl sein!'); +end; + +if (C ~= fix(C)), + error('Fehler: C muss eine ganze Zahl sein!'); +end; + +if (P ~= fix(P)), + error('Fehler: P muss eine ganze Zahl sein!'); +end; + +if (S ~= fix(S)), + error('Fehler: S muss eine ganze Zahl sein!'); +end; + +if (L ~= fix(L)), + error('Fehler: L muss eine ganze Zahl sein!'); +end; + +% ------------------------------------------------------------------------- +% Automatik-Auswahl +% ------------------------------------------------------------------------- +% Falls P=0 => Automatische Wahl von P +% Anahme: L+N/P = 2N/P +if (P==0), + P = N/L; + if (rem(N,L) ~= 0) + error('Fehler: Quotient aus N und L muss eine ganze Zahl sein! (N = P*S*L)'); + end; +end; + +% Falls C=0 => Automatische Wahl von C=2^r +% Anahme: C = L + N/P +if (C==0) + C = 2^ceil(log2(L+N/P-1)); +end; + +% ------------------------------------------------------------------------- +% FLMS Parameter +% ------------------------------------------------------------------------- +gamma = 1.0/C; % Sicherheitskonstante vermeidet Division durch Null +lambda = 0.4; % Vergessensfaktor + +% ------------------------------------------------------------------------- +% Initialisierungen +% ------------------------------------------------------------------------- +Lx = length(x); % Länge des Eingangsvektors +Np = S*L; % Teilfilter hat die Länge Np = S*L = N/P +K = Lx/L; % Anzahl Verarbeitungsblöcke + +e=zeros(Lx,1); % Fehlervektor +dd=zeros(Lx,1); % Filterausgang, Schaetzung von d +wp = zeros(Np,P); % P Teilfilter im Zeitbereich +PX = zeros(C,1); % Schaetzung des Leistungsdichtespektrums + +WSp = zeros(C,P); % P Teilfilter im Frequenzbereich +X = zeros(C,S*P); % Eingangsvektor im Frequenzbereich +mu_px = zeros(C,P*S); % Variable Schrittweite +dd_ = zeros(C,1); % Filterausgang im Zeitbereich +DD = zeros(C,1); % Filterausgang im Frequenzbereich +xzp = zeros(Lx+S*L,1); % Eingangsvektor mit S*L führenden Nullen + + +% ------------------------------------------------------------------------- +% Tabellen zur effizienten MatLab-Implementierung +% ------------------------------------------------------------------------- +pidTbl = [2:P*S]; % Index auf P*S vergangene Eingangsvektoren +pidTbl(P*S) = 1; % Reihenfolge: x(k-1),x(k-2),..,x(k-P*S-1),x(k) + +pStbl = [0:P*S-1]; % Tabelle zur korrekten Adressierung der Eingangsvektoren... +pStbl(1) = P*S; % ...bei der Filterung und Update der Koeffizienten. + +pjiTbl = [2:P]; % Index auf Teilfilter, welches bei der effizienten... +pjiTbl(P) = 1; % ...Projektion bearbeitet wird + +xzp(C-L+1:Lx+C-L) = x; % Einfügen von Np führenden Nullen an den Eingangsvektor + +pid = 1; % Aktueller Index zeigt auf Eingangsvektor x(0) +pji = 1; % Aktueller Index zeigt auf Teilfilter w_0 + + +% ------------------------------------------------------------------------- +% Anfangswerte der Filterkoeffienten partitionieren +% und in den Frequenzbereich transformieren wi_p(k) => WS_p(k) +% ------------------------------------------------------------------------- +for p=1:P, + SLp = (p-1)*N/P; + wp(1:Np,p) = wi(SLp+1:SLp+N/P); + WSp(1:C,p) = 1/C*fft(wp(1:Np,p),C); +end; + +% ------------------------------------------------------------------------- +% Algorithmus Start +% ------------------------------------------------------------------------- +flops(0); + +for k=1:K, + + kL = (k-1)*L; + + % Transformation des aktuellen Eingangsvektors in den Frequenzbereich + X(1:C,pid) = 1/C*fft(xzp(kL+1:kL+C),C); + + % Partitioned FFT-Faltung + % Faltung des ersten Teilfilters mit aktuellem Eingangsvektor + pS = pid; + DD = WSp(1:C,1).*X(1:C,pid) *C; + + % Faltung der restlichen Teilfilters mit vergangenen Eingangsvektoren + for p=2:P, + % Ermittlung des korrekten Index 'pS' aus Tabelle + for i = 1:S, + pS = pStbl(pS); + end; + % Summation der einzelnen Filterausgänge + DD = DD + WSp(1:C,p).*X(1:C,pS) *C; + end; + + + % Transformation des Filterausgangs in den Zeitbereich + dd_ = real(ifft(DD,C))*C; + + % Abspeichern der letzten L Werte des Filterausgangs + dd(kL+1:kL+L) = dd_(C-L+1:C); + + % Berechnung des Fehlers e(k) + e(kL+1:kL+L) = d(kL+1:kL+L) - dd(kL+1:kL+L); + + % Transformation des Fehlers in den Frequenzbereich + E = 1/C*fft([zeros(C-L,1); e(kL+1:kL+L).*mu(kL+1:kL+L)],C); + + % Schätzung der mittleren Eingangsleistung Px(k) aus x(k) + PX = abs(lambda*conj(X(1:C,pid)).*X(1:C,pid) *C + (1-lambda)*PX); + + + % Berechnung der variablen Schrittweite mu(k) + % mu(k) wird niemals grösser als eins + mu_px(1:C,pid) = (gamma/P)./(PX+gamma); + + % Filterkoeffizenten-Update + pS = pid; + for p=1:P, + % Aktualisierung des p-ten Teilfilters + WSp(1:C,p) = WSp(1:C,p) + mu_px(1:C,pS) .* conj(X(1:C,pS)) .* E *C; + % Ermittlung des korrekten Index 'pS' aus Tabelle + for i = 1:S, + pS = pStbl(pS); + end; + end; + +% Aufwändig: Projektion jedes p-ten Teilfilters wp pro k-ter Iteration +% for p=1:P, +% wp(1:C,p) = real(dspifft(WSp(1:C,p),C)); +% WSp(1:C,p) = dspfft(wp(1:Np,p),C); +% end; + +% Effizient: Projektion eines pji-ten Teilfilters alternierend pro k-ter Iteration + wp(1:C,pji) = C*real(ifft(WSp(1:C,pji),C)); % pji = 0,1,..,P-1,0,1,..P-1,.. + WSp(1:C,pji) = 1/C*fft(wp(1:Np,pji),C); + + pji = pjiTbl(pji); % Nächster Index pji auf Eingangsvektor aus Tabelle + pid = pidTbl(pid); % Nächster Index pid auf Teilfilter aus Tabelle + +end; + +% Komposition des Gesamtfilters aus den P Teilfiltern +for p=1:P, + wf((p-1)*Np+1:p*Np) = wp(1:Np,p); +end; + +fpo=flops; + +% ------------------------------------------------------------------------- +% Ende PFLMS.M +% ------------------------------------------------------------------------- diff --git a/common/pflms.old.m b/common/pflms.old.m new file mode 100755 index 0000000..d823168 --- /dev/null +++ b/common/pflms.old.m @@ -0,0 +1,211 @@ +% Partitioned Frequency LMS Adaptive Filter -PFLMS- +% ------------------------------------------------- +% Aufruf: +% Function [ys,ee,wf,fpo]=pflms(x,d,alpha,N,C,P,S,L,wi) +% +% Parameter: +% x : Eingangsvektor +% d : Referenzvektor +% alpha: konstante Schrittweite 0 < alpha <= 1.0 +% N : Anzahl Filterkoeffizienten (Hinweis: N = P*S*L) +% C : FFT-Länge, idealerweise eine Zweier-Potenz (0=Auto) +% P : Anzahl der Filterpartitionen +% S : Anzahl der Filtersegmente +% L : Länge eines Verarbeitungsblocks (idealerweise C=L+N/P=2N/P) +% wi : Startwerte der Filterkoeffizienten im Zeitbereich +% +% Rückgabewerte: +% ee : Adaptitionsfehler +% ys : Filterausgang +% wf : Filterkoeffizienten am Ende der Adaption im Zeitbereich +% fpo : Anzahl benötigten der Floating-Point Operationen + +% ------------------------------------------------------------------------- +% Datum : 20.7.2001 +% Autor : Jens Ahrensfeld +% Thema : Diplomarbeit +% Datei : pflms.m +% Benötigte Dateien: dspfft.m, dspifft.m +% +% ------------------------------------------------------------------------- +function [ys,ee,wf,fpo]=pflms(x,d,alpha,N,C,P,S,L,wi) + +% ------------------------------------------------------------------------- +% Überprüfen der Parameter +% ------------------------------------------------------------------------- +if (N ~= fix(N)), + error('Fehler: N muss eine ganze Zahl sein!'); +end; + +if (C ~= fix(C)), + error('Fehler: C muss eine ganze Zahl sein!'); +end; + +if (P ~= fix(P)), + error('Fehler: P muss eine ganze Zahl sein!'); +end; + +if (S ~= fix(S)), + error('Fehler: S muss eine ganze Zahl sein!'); +end; + +if (L ~= fix(L)), + error('Fehler: L muss eine ganze Zahl sein!'); +end; + +% ------------------------------------------------------------------------- +% Automatik-Auswahl +% ------------------------------------------------------------------------- +% Falls P=0 => Automatische Wahl von P +% Anahme: L+N/P = 2N/P +if (P==0), + P = N/L; + if (rem(N,L) ~= 0) + error('Fehler: Quotient aus N und L muss eine ganze Zahl sein! (N = P*S*L)'); + end; +end; + +% Falls C=0 => Automatische Wahl von C=2^r +% Anahme: C = L + N/P +if (C==0) + C = 2^ceil(log2(L+N/P-1)); +end; + +% ------------------------------------------------------------------------- +% FLMS Parameter +% ------------------------------------------------------------------------- +gamma = 1.0/C; % Sicherheitskonstante vermeidet Division durch Null +lambda = 0.6; % Vergessensfaktor + +% ------------------------------------------------------------------------- +% Initialisierungen +% ------------------------------------------------------------------------- +Lx = length(x); % Länge des Eingangsvektors +Np = S*L; % Teilfilter hat die Länge Np = S*L = N/P +K = Lx/L; % Anzahl Verarbeitungsblöcke + +ee=zeros(Lx,1); % Fehlervektor +ys=zeros(Lx,1); % Filterausgang, Schaetzung von d +wp = zeros(Np,P); % P Teilfilter im Zeitbereich +PX = zeros(C,1); % Schaetzung des Leistungsdichtespektrums + +WSp = zeros(C,P); % P Teilfilter im Frequenzbereich +X = zeros(C,S*P); % Eingangsvektor im Frequenzbereich +mu = zeros(C,P*S); % Variable Schrittweite +yy = zeros(C,1); % Filterausgang im Zeitbereich +YS = zeros(C,1); % Filterausgang im Frequenzbereich +xzp = zeros(Lx+S*L,1); % Eingangsvektor mit S*L führenden Nullen + + +% ------------------------------------------------------------------------- +% Tabellen zur effizienten MatLab-Implementierung +% ------------------------------------------------------------------------- +pidTbl = [2:P*S]; % Index auf P*S vergangene Eingangsvektoren +pidTbl(P*S) = 1; % Reihenfolge: x(k-1),x(k-2),..,x(k-P*S-1),x(k) + +pStbl = [0:P*S-1]; % Tabelle zur korrekten Adressierung der Eingangsvektoren... +pStbl(1) = P*S; % ...bei der Filterung und Update der Koeffizienten. + +pjiTbl = [2:P]; % Index auf Teilfilter, welches bei der effizienten... +pjiTbl(P) = 1; % ...Projektion bearbeitet wird + +xzp(C-L+1:Lx+C-L) = x; % Einfügen von Np führenden Nullen an den Eingangsvektor + +pid = 1; % Aktueller Index zeigt auf Eingangsvektor x(0) +pji = 1; % Aktueller Index zeigt auf Teilfilter w_0 + + +% ------------------------------------------------------------------------- +% Anfangswerte der Filterkoeffienten partitionieren +% und in den Frequenzbereich transformieren wi_p(k) => WS_p(k) +% ------------------------------------------------------------------------- +for p=1:P, + SLp = (p-1)*N/P; + wp(1:Np,p) = wi(SLp+1:SLp+N/P); + WSp(1:C,p) = 1/C*fft(wp(1:Np,p),C); +end; + +% ------------------------------------------------------------------------- +% Algorithmus Start +% ------------------------------------------------------------------------- +flops(0); + +for k=1:K, + + kL = (k-1)*L; + + % Transformation des aktuellen Eingangsvektors in den Frequenzbereich + X(1:C,pid) = 1/C*fft(xzp(kL+1:kL+C),C); + + % Partitioned FFT-Faltung + % Faltung des ersten Teilfilters mit aktuellem Eingangsvektor + pS = pid; + YS = WSp(1:C,1).*X(1:C,pid) *C; + + % Faltung der restlichen Teilfilters mit vergangenen Eingangsvektoren + for p=2:P, + % Ermittlung des korrekten Index 'pS' aus Tabelle + for i = 1:S, + pS = pStbl(pS); + end; + % Summation der einzelnen Filterausgänge + YS = YS + WSp(1:C,p).*X(1:C,pS) *C; + end; + + + % Transformation des Filterausgangs in den Zeitbereich + yy = real(ifft(YS,C))*C; + + % Abspeichern der letzten L Werte des Filterausgangs + ys(kL+1:kL+L) = yy(C-L+1:C); + + % Berechnung des Fehlers e(k) + ee(kL+1:kL+L) = d(kL+1:kL+L) - ys(kL+1:kL+L); + + % Transformation des Fehlers in den Frequenzbereich + E = 1/C*fft([zeros(C-L,1); ee(kL+1:kL+L)],C); + + % Schätzung der mittleren Eingangsleistung Px(k) aus x(k) + PX = abs((1-lambda)*conj(X(1:C,pid)).*X(1:C,pid) *C + lambda*PX); + + + % Berechnung der variablen Schrittweite mu(k) + % mu(k) wird niemals grösser als eins + mu(1:C,pid) = (alpha*gamma/P) ./(PX+gamma); + + % Filterkoeffizenten-Update + pS = pid; + for p=1:P, + % Aktualisierung des p-ten Teilfilters + WSp(1:C,p) = WSp(1:C,p) + mu(1:C,pS) .* conj(X(1:C,pS)) .* E *C; + % Ermittlung des korrekten Index 'pS' aus Tabelle + for i = 1:S, + pS = pStbl(pS); + end; + end; + +% Aufwändig: Projektion jeden p-ten Teilfilters wp pro k-ter Iteration +% for p=1:P, +% wp(1:C,p) = real(dspifft(WSp(1:C,p),C)); +% WSp(1:C,p) = dspfft(wp(1:Np,p),C); +% end; + +% Effizient: Projektion eines Teilfilters alternierend pro k-ter Iteration + wp(1:C,pji) = C*real(ifft(WSp(1:C,pji),C)); + WSp(1:C,pji) = 1/C*fft(wp(1:Np,pji),C); + + pji = pjiTbl(pji); % Nächster Index auf Eingangsvektor aus Tabelle + pid = pidTbl(pid); % Nächster Index auf Teilfilter aus Tabelle + +end; + +% Komposition des Gesamtfilters aus den P Teilfiltern +for p=1:P, + wf((p-1)*Np+1:p*Np) = wp(1:Np,p); +end; + +fpo=flops; + +% ------------------------------------------------------------------------- +% Ende PFLMS.M +% ------------------------------------------------------------------------- diff --git a/common/pflms2.m b/common/pflms2.m new file mode 100755 index 0000000..7818455 --- /dev/null +++ b/common/pflms2.m @@ -0,0 +1,97 @@ +% pflms2(mu,N,C,P,S,L,x,d,w_start) +function [ee,w,dd,PX,u]=pflms2(mu,N,C,P,S,L,x,d,w_start) + +Lx = length(x); +Np = S*L; +K = Lx/L; + +% FLMS Parameter +gamma = 0.6; % Vergessensfaktor +alpha = 1.0 % 0 < alpha < 1 +beta = mu/C + +% Initialisierungen +ee=zeros(Lx,1); % Fehlervektor +dd=zeros(Lx,1); % Filterausgang, Schaetzung von d +wp = zeros(C,P); % Gewichtsvektor im Frequenzbereich +w = zeros(P*S*L,1); % Gewichtsvektor im Zeitbereich +PX = zeros(C,1); % Schaetzung des Leistungsdichtespektrums + +WSp = zeros(C,P); +X = zeros(C,S*P); +u = zeros(C,P*S); +ys = zeros(C,1); +YS = zeros(C,1); +xzp = zeros(Lx+S*L,1); +pidTbl = [2:P*S]; +pidTbl(P*S) = 1; +pStbl = [0:P*S-1]; +pStbl(1) = P*S; +pjiTbl = [2:P]; +pjiTbl(P) = 1; + +xzp(C-L+1:Lx+C-L) = x; + +pid = 1; +pji = 1; + +flops(0); + +for k=1:K, + + kL = (k-1)*L; + + X(1:C,pid) = dspfft(xzp(kL+1:kL+C),C); + + pS = pid; + YS = WSp(1:C,1).*X(1:C,pid) *C; + + for p=2:P, + for i = 1:S + pS = pStbl(pS); + end; + YS = YS + WSp(1:C,p).*X(1:C,pS) *C; + end; + + ys = real(dspifft(YS,C)); + dd(kL+1:kL+L) = ys(C-L+1:C); + + ee(kL+1:kL+L) = d(kL+1:kL+L) - ys(C-L+1:C); + E = dspfft([zeros(C-L,1); ee(kL+1:kL+L)],C); + + PX = abs((1-gamma)*conj(X(1:C,pid)).*X(1:C,pid) *C + gamma*PX); + + % umax = alpha + u(1:C,pid) = (alpha*beta) ./(PX+beta); + + % Update + pS = pid; + for p=1:P, + WSp(1:C,p) = WSp(1:C,p) + u(1:C,pS) .* conj(X(1:C,pS)) .* E *C; + for i = 1:S + pS = pStbl(pS); + end; + end; + +% Teuer: Projektion jeden p-ten Teilfilters wp pro k-ter Iteration +% for p=1:P, +% wp(1:C,p) = real(dspifft(WSp(1:C,p),C)); +% WSp(1:C,p) = dspfft(wp(1:Np,p),C); +% end; + +% Billig: Projektion eines Teilfilters alternierend pro k-ter Iteration + wp(1:C,pji) = real(dspifft(WSp(1:C,pji),C)); + WSp(1:C,pji) = dspfft(wp(1:Np,pji),C); + + pji = pjiTbl(pji); + pid = pidTbl(pid); + +end; + +flops + +% Return estimated filter weights +for p=1:P, + w((p-1)*Np+1:p*Np) = wp(1:Np,p); +end; + diff --git a/common/pflms_d.m b/common/pflms_d.m new file mode 100755 index 0000000..447a828 --- /dev/null +++ b/common/pflms_d.m @@ -0,0 +1,219 @@ +% Partitioned Frequency LMS Adaptive Filter -PFLMS- (Debug Version) +% ----------------------------------------------------------------- +% Aufruf: +% Function [ys,ee,wf,fpo]=pflms_d(x,d,alpha,N,C,P,S,L,wi) +% +% Parameter: +% x : Eingangsvektor +% d : Referenzvektor +% alpha: konstante Schrittweite 0 < alpha <= 1.0 +% N : Anzahl Filterkoeffizienten (Hinweis N = P*S*L) +% C : FFT-Länge (idealerweise eine Zweier-Potenz) +% P : Anzahl der Filterpartitionen +% S : Anzahl der Filtersegmente +% L : Länge eines Verarbeitungsblocks (idealerweise C=L+N/P=2N/P) +% wi : Startwerte der Filterkoeffizienten im Zeitbereich +% +% Rückgabewerte: +% ee : Adaptitionsfehler +% ys : Filterausgang +% wf : Filterkoeffizienten am Ende der Adaption im Zeitbereich +% fpo : Anzahl benötigten der Floating-Point Operationen +% +% Zu Debug-Zwecken werden folgende Variablen als Binär-Dateien gespeichert: +% wf.dat : Filterkoeffizienten am Ende der Adaption (Länge N) + +% ------------------------------------------------------------------------- +% Datum : 20.7.2001 +% Autor : Jens Ahrensfeld +% Thema : Diplomarbeit +% Datei : pflms_d.m +% Benötigte Dateien: dspfft.m, dspifft.m +% +% ------------------------------------------------------------------------- +function [ys,ee,wf,fpo]=pflms_d(x,d,alpha,N,C,P,S,L,wi) + +% ------------------------------------------------------------------------- +% Überprüfen der Parameter +% ------------------------------------------------------------------------- +if (N ~= fix(N)), + error('Fehler: N muss eine ganze Zahl sein!'); +end; + +if (C ~= fix(C)), + error('Fehler: C muss eine ganze Zahl sein!'); +end; + +if (P ~= fix(P)), + error('Fehler: P muss eine ganze Zahl sein!'); +end; + +if (S ~= fix(S)), + error('Fehler: S muss eine ganze Zahl sein!'); +end; + +if (L ~= fix(L)), + error('Fehler: L muss eine ganze Zahl sein!'); +end; + +% ------------------------------------------------------------------------- +% Automatik-Auswahl +% ------------------------------------------------------------------------- +% Falls P=0 => Automatische Wahl von P +% Anahme: L+N/P = 2N/P +if (P==0), + P = N/L + if (rem(N,L) ~= 0) + error('Fehler: Quotient aus N und L muss eine ganze Zahl sein! (N = P*S*L)'); + end; +end; + +% Falls C=0 => Automatische Wahl von C=2^r +% Anahme: C = L + N/P +if (C==0) + C = 2^ceil(log2(L+N/P-1)) +end; + +% ------------------------------------------------------------------------- +% FLMS Parameter +% ------------------------------------------------------------------------- +gamma = 1.0/C; % Sicherheitskonstante vermeidet Division durch Null +lambda = 0.6; % Vergessensfaktor + +% ------------------------------------------------------------------------- +% Initialisierungen +% ------------------------------------------------------------------------- +Lx = length(x); % Länge des Eingangsvektors +Np = S*L; % Teilfilter hat die Länge Np = S*L = N/P +K = Lx/L; % Anzahl Verarbeitungsblöcke + +ee=zeros(Lx,1); % Fehlervektor +ys=zeros(Lx,1); % Filterausgang, Schaetzung von d +wp = zeros(C,P); % P Teilfilter im Zeitbereich +wf = zeros(P*S*L,1); % Gesamtfilter im Zeitbereich +PX = zeros(C,1); % Schaetzung des Leistungsdichtespektrums + +WSp = zeros(C,P); % P Teilfilter im Frequenzbereich +X = zeros(C,S*P); % Eingangsvektor im Frequenzbereich +mu = zeros(C,P*S); % Variable Schrittweite +yy = zeros(C,1); % Filterausgang im Zeitbereich +YS = zeros(C,1); % Filterausgang im Frequenzbereich +xzp = zeros(Lx+S*L,1); % Eingangsvektor mit S*L führenden Nullen + +% ------------------------------------------------------------------------- +% Tabellen zur effizienten MatLab-Implementierung +% ------------------------------------------------------------------------- +pidTbl = [2:P*S]; % Index auf P*S vergangene Eingangsvektoren +pidTbl(P*S) = 1; % Reihenfolge: x(k-1),x(k-2),..,x(k-P*S-1),x(k) + +pStbl = [0:P*S-1]; % Tabelle zur korrekten Adressierung der Eingangsvektoren... +pStbl(1) = P*S; % ...bei der Filterung und Update der Koeffizienten. + +pjiTbl = [2:P]; % Index auf Teilfilter, welches bei der effizienten... +pjiTbl(P) = 1; % ...Projektion bearbeitet wird + +xzp(C-L+1:Lx+C-L) = x; % Einfügen von Np führenden Nullen an den Eingangsvektor + +% Start +pid = 1; % Aktueller Index zeigt auf Eingangsvektor x(0) +pji = 1; % Aktueller Index zeigt auf Teilfilter w_0 + +% ------------------------------------------------------------------------- +% Anfangswerte der Filterkoeffienten partitionieren +% und in den Frequenzbereich transformieren wi_p(k) => WS_p(k) +% ------------------------------------------------------------------------- +%for p=1:P, +% SLp = (p-1)*N/P; +% wp(p,1:N/P) = wi(SLp+1:SLp+N/P); +% WSp(p,1:C) = dspfft(wp(p,1:C),C); +%end; + + +% ------------------------------------------------------------------------- +% Algorithmus Start +% ------------------------------------------------------------------------- +flops(0); + +for k=1:K, + + kL = (k-1)*L; + + % Transformation des aktuellen Eingangsvektors in den Frequenzbereich + X(1:C,pid) = dspfft(xzp(kL+1:kL+C),C); + + % Partitioned FFT-Faltung + % Faltung des ersten Teilfilters mit aktuellem Eingangsvektor + pS = pid; + YS = WSp(1:C,1).*X(1:C,pid) *C; + + % Faltung der restlichen Teilfilters mit vergangenen Eingangsvektoren + for p=2:P, + % Ermittlung des korrekten Index 'pS' aus Tabelle + for i = 1:S, + pS = pStbl(pS); + end; + % Summation der einzelnen Filterausgänge + YS = YS + WSp(1:C,p).*X(1:C,pS) *C; + end; + + + % Transformation des Filterausgangs in den Zeitbereich + yy = real(dspifft(YS,C)); + + % Abspeichern der letzten L Werte des Filterausgangs + ys(kL+1:kL+L) = yy(C-L+1:C); + + % Berechnung des Fehlers e(k) + ee(kL+1:kL+L) = d(kL+1:kL+L) - ys(kL+1:kL+L); + + % Transformation des Fehlers in den Frequenzbereich + E = dspfft([zeros(C-L,1); ee(kL+1:kL+L)],C); + + % Schätzung der mittleren Eingangsleistung Px(k) aus x(k) + PX = abs((1-lambda)*conj(X(1:C,pid)).*X(1:C,pid) *C + lambda*PX); + + + % Berechnung der variablen Schrittweite mu(k) + % mu(k) wird niemals grösser als eins + mu(1:C,pid) = (alpha*gamma/P) ./(PX+gamma); + + % Filterkoeffizenten-Update + pS = pid; + for p=1:P, + % Aktualisierung des p-ten Teilfilters + WSp(1:C,p) = WSp(1:C,p) + mu(1:C,pS) .* conj(X(1:C,pS)) .* E *C; + % Ermittlung des korrekten Index 'pS' aus Tabelle + for i = 1:S, + pS = pStbl(pS); + end; + end; + +% Aufwändig: Projektion jeden p-ten Teilfilters wp pro k-ter Iteration +% for p=1:P, +% wp(1:C,p) = real(dspifft(WSp(1:C,p),C)); +% WSp(1:C,p) = dspfft(wp(1:Np,p),C); +% end; + +% Effizient: Projektion eines Teilfilters alternierend pro k-ter Iteration + wp(1:C,pji) = real(dspifft(WSp(1:C,pji),C)); + WSp(1:C,pji) = dspfft(wp(1:Np,pji),C); + + pji = pjiTbl(pji); % Nächster Index auf Eingangsvektor aus Tabelle + pid = pidTbl(pid); % Nächster Index auf Teilfilter aus Tabelle + +end; + +% Komposition des Gesamtfilters aus den P Teilfiltern +for p=1:P, + wf((p-1)*Np+1:p*Np) = wp(1:Np,p); +end; + +fpo=flops; + +% Speichern der Koeffizienten in Datei (Genauigkeit Float32) +fid = fopen('wf.dat','wb'); +fwrite(fid,wf,'float32'); +fclose(fid); + +% ------------------------------------------------------------------------- +% Ende PFLMS_D.M diff --git a/common/pflmss.m b/common/pflmss.m new file mode 100755 index 0000000..7bb7e44 --- /dev/null +++ b/common/pflmss.m @@ -0,0 +1,96 @@ +% pflmss(mu,N,C,P,S,L,x,d,w_start) +function [ee,w,dd,calc,PX]=pflmss(mu,N,C,P,S,L,x,d,w_start) + +Lx = length(x); +Np = S*L; +K = Lx/L; + +% Initialisierungen + +%% FLMS Parameter +gamma = 0.6; % Vergessensfaktor + +ee=zeros(Lx,1); % Fehlervektor +dd=zeros(Lx,1); % Filterausgang, Schaetzung von d +wp = zeros(C,P); % Gewichtsvektor im Frequenzbereich +w = zeros(P*S*L,1); % Gewichtsvektor im Zeitbereich +PX = ones(C,1); % Schaetzung des Leistungsdichtespektrums + +WSp = zeros(C,P); +X = zeros(C,S*P); +U = zeros(C,S*P); +ys = zeros(C,1); +YS = zeros(C,1); +xzp = zeros(Lx+S*L,1); +pidTbl = [2:P*S]; +pidTbl(P*S) = 1; +pStbl = [0:P*S-1]; +pStbl(1) = P*S; +pjiTbl = [2:P]; +pjiTbl(P) = 1; + +xzp(C-L+1:Lx+C-L) = x; + +tic; % Stopuhr laeuft + +pid = 1; +pji = 1; + +flops(0); + +for k=1:K, + + kL = (k-1)*L; + + X(1:C,pid) = fft(xzp(kL+1:kL+C),C); + + pS = pid; + YS = WSp(1:C,1).*X(1:C,pid); + for p=2:P, + for i = 1:S + pS = pStbl(pS); + end; + YS = YS + WSp(1:C,p).*X(1:C,pS); + end; + + ys = real(ifft(YS,C)); + dd(kL+1:kL+L) = ys(C-L+1:C); + + ee(kL+1:kL+L) = d(kL+1:kL+L) - ys(C-L+1:C); + E = fft([zeros(C-L,1); ee(kL+1:kL+L)]); + PX = abs((1-gamma)*conj(X(1:C,pid)).*X(1:C,pid) + gamma*PX); + + U(1:C,pid) = mu ./ (PX+0.001); + + % Update + pS = pid; + for p=1:P, + WSp(1:C,p) = WSp(1:C,p) + U(1:C,pS) .* conj(X(1:C,pS)) .* E; + for i = 1:S + pS = pStbl(pS); + end; + end; + +% Teuer: Projektion jeden p-ten Teilfilters wp pro k-ter Iteration +% for p=1:P, +% wp(1:C,p) = real(ifft(WSp(1:C,p))); +% WSp(1:C,p) = fft(wp(1:Np,p),C); +% end; + +% Billig: Projektion eines Teilfilters alternierend pro k-ter Iteration + wp(1:C,pji) = real(ifft(WSp(1:C,pji))); + WSp(1:C,pji) = fft(wp(1:Np,pji),C); + + pji = pjiTbl(pji); + pid = pidTbl(pid); + +end; + +flops + +% Return estimated filter weights +for p=1:P, + w((p-1)*Np+1:p*Np) = wp(1:Np,p); +end; + +calc = toc; \ No newline at end of file diff --git a/common/phnl.m b/common/phnl.m new file mode 100755 index 0000000..aea39d9 --- /dev/null +++ b/common/phnl.m @@ -0,0 +1,15 @@ +% [nl,nlf] = phnl(e,es,el,xs,xl,al,pl,nli) +function [nl,nlf] = phnl(e,es,el,xs,xl,al,pl,nli) + +N = lge(e); +e = abs(e); +nl(1) = nli; + +for n = 2:N, + if((xs(n) <= pl*xl(n-1)) & (es(n) <= pl*el(n))) + nl(n) = al*e(n) + (1-al)*nl(n-1); + else + nl(n) = nl(n-1); + end; +end; +nlf = nl(N); diff --git a/common/phxl.m b/common/phxl.m new file mode 100755 index 0000000..8e36233 --- /dev/null +++ b/common/phxl.m @@ -0,0 +1,15 @@ +% [xl,xlf] = phxl(x,xs,al,pl,xli) +function [xl,xlf] = phxl(x,xs,al,pl,xli) + +N = lge(x); +x = abs(x); +xl(1) = xli; + +for n = 2:N, + if(xs(n) > pl*xl(n-1)) + xl(n) = xl(n-1); + else + xl(n) = al*x(n) + (1-al)*xl(n-1); + end; +end; +xlf = xl(N); diff --git a/common/phxs.m b/common/phxs.m new file mode 100755 index 0000000..2f36284 --- /dev/null +++ b/common/phxs.m @@ -0,0 +1,15 @@ +% [xs,xsf] = phxs(x,ar,af,xsi) +function [xs,xsf] = phxs(x,ar,af,xsi) + +N = lge(x); +x = abs(x); +xs(1) = xsi; + +for n = 2:N, + if(x(n) > xs(n-1)) + xs(n) = ar*x(n) + (1-ar)*xs(n-1); + else + xs(n) = af*x(n) + (1-af)*xs(n-1); + end; +end; +xsf = xs(N); diff --git a/common/pltimp.m b/common/pltimp.m new file mode 100755 index 0000000..ef6ccad --- /dev/null +++ b/common/pltimp.m @@ -0,0 +1,66 @@ +% pltimp(x,Fa,Nw,NP) +function y = pltimp(x,Fa,Nw,NP) + +xn = x(1:NP)/max(x(1:NP)); +Px = sm2(xn(1:NP).^2,Nw,0); +Pxn = Px/max(Px(1:NP)); +PxndB = 10*log(Pxn(1:NP)); + +matchval = -120; +for n = 1:NP, + if ((PxndB(n) <= -60) & (PxndB(n) > matchval)) + matchn = n; + matchval = PxndB(n); + end; + if (PxndB(n) > -60) + matchval = -120; + end; +end; + +matchtime = matchn*1000/Fa +matchval + +figure; + plot(1000*(1:NP)/Fa,x(1:NP)); + grid; + PlotTitle = sprintf('h(t) Fa=%d Hz',Fa); + title(PlotTitle); + xlabel('t [ms]'); + ylabel('h'); + +figure; + plot(1000*(1:NP)/Fa,Px(1:NP)); + grid; + PlotTitle = sprintf('Leistung Ph(t), Windowsize=%d Samples, Time=%g ms',Nw,1000*Nw/Fa); + title(PlotTitle); + xlabel('t [ms]'); + ylabel('eh'); + +figure; + plot(1000*(1:NP)/Fa,PxndB(1:NP)); + set(gca,'YTick',-20*(0:6)); + xt = get(gca,'XTick'); + rval = 100; + di = 1; + for i = 1:length(xt)-1, + lv = round(matchtime) - xt(i); + if (lv <= 0) + di = i; + if (i > 1) + uv = round(matchtime) - xt(i-1); + if (abs(uv) < abs(lv)) + di = i - 1; + end; + break; + end; + end; + end; + di + xt(di) = matchn*1000/Fa; + set(gca,'XTick',xt); + grid; + PlotTitle = sprintf('Leistung Ph(t) 60dB-decay time=%g ms',matchtime); + title(PlotTitle); + xlabel('t [ms]'); + ylabel('Ph [dB]'); + \ No newline at end of file diff --git a/common/qamTable.m b/common/qamTable.m new file mode 100755 index 0000000..cb47610 --- /dev/null +++ b/common/qamTable.m @@ -0,0 +1,17 @@ +% qam(n) +% Qam +% + +function qamTbl=qamTable(n) +startIQ = 1.0/sqrt(2.0) +stepIQ = 2*startIQ/(sqrt(n)-1) +qamTbl = zeros(n,1); + +s=1; +for I=-startIQ:stepIQ:startIQ, + for Q=-startIQ:stepIQ:startIQ, + qamTbl(s) = I + i*Q; + s = s + 1; + end; +end; + diff --git a/common/radio/fir_rc.bak.m b/common/radio/fir_rc.bak.m new file mode 100755 index 0000000..bf18076 --- /dev/null +++ b/common/radio/fir_rc.bak.m @@ -0,0 +1,24 @@ +% function [b,k] = fir_rc(fa, T, N, alpha) +% Raised Cosine FIR-Filter + +function [b,k] = fir_rc(fa, T, alpha, N) + +n = -(N-1)/2:1:(N-1)/2; +t = n/fa; + +nom = sinc(t/T).*cos(pi*alpha*t/T); +denom = (1 - (2*alpha*t/T).^2); + +% Find zero elements +ze = find(abs(alpha*t/T) == 0.5); + +nom(ze) = pi/4 * sinc(t(ze)./T); +denom(ze) = 1; + +b = (nom./denom).*blackman(N)'; +k = 1 /sum(b); + + + + + diff --git a/common/radio/fir_rc.m b/common/radio/fir_rc.m new file mode 100755 index 0000000..4d8881f --- /dev/null +++ b/common/radio/fir_rc.m @@ -0,0 +1,25 @@ +% [h] =fir_rc(Fa,Tsym,Alpha,N) +function [b,k,delay] =fir_rc(fa,Tsym,Alpha,N) + +if rem(N,2), + delay = (N-1)/2; +else + delay = N/2; +end + +n = ((0:N-1)-delay)/fa; + +ind1 = find(abs(abs(4.*Alpha.*n./Tsym) - 1.0) > sqrt(eps)); +if ~isempty(ind1), + nind = n(ind1); + b(ind1) = sinc(2.*nind./Tsym)./fa ... + .* cos(2.*pi.*Alpha.*nind./Tsym) ... + ./ (1.0 - (4.*Alpha.*nind./Tsym).^2); +end + +ind = 1:length(n); +ind(ind1) = []; +b(ind) = Alpha ./ (2.*fa) .* sin(pi ./ (2.*Alpha)); + +b = 2.*b./Tsym; +k=0.5*fa*Tsym; diff --git a/common/radio/fir_rc_bak.m b/common/radio/fir_rc_bak.m new file mode 100755 index 0000000..bf18076 --- /dev/null +++ b/common/radio/fir_rc_bak.m @@ -0,0 +1,24 @@ +% function [b,k] = fir_rc(fa, T, N, alpha) +% Raised Cosine FIR-Filter + +function [b,k] = fir_rc(fa, T, alpha, N) + +n = -(N-1)/2:1:(N-1)/2; +t = n/fa; + +nom = sinc(t/T).*cos(pi*alpha*t/T); +denom = (1 - (2*alpha*t/T).^2); + +% Find zero elements +ze = find(abs(alpha*t/T) == 0.5); + +nom(ze) = pi/4 * sinc(t(ze)./T); +denom(ze) = 1; + +b = (nom./denom).*blackman(N)'; +k = 1 /sum(b); + + + + + diff --git a/common/radio/fir_rrc.m b/common/radio/fir_rrc.m new file mode 100755 index 0000000..0209147 --- /dev/null +++ b/common/radio/fir_rrc.m @@ -0,0 +1,35 @@ +% [h] =fir_rrc(fa,Tsym,Alpha,N) +function [b,k,delay] =fir_rrc(fa,Tsym,Alpha,N) + +if rem(N,2), + delay = (N-1)/2; +else + delay = N/2; +end + +n = ((0:N-1)-delay)/fa; + +ind1 = find(n == 0); +if ~isempty(ind1), + b(ind1) = - sqrt(2./Tsym) ./ (pi.*fa) .* (pi.*(Alpha-1) - 4.*Alpha ); +end + +ind2 = find(abs(abs(8.*Alpha.*n./Tsym) - 1.0) < sqrt(eps)); +if ~isempty(ind2), + b(ind2) = sqrt(2./Tsym) ./ (2.*pi.*fa) ... + * ( pi.*(Alpha+1) .* sin(pi.*(Alpha+1)./(4.*Alpha)) ... + - 4.*Alpha .* sin(pi.*(Alpha-1)./(4.*Alpha)) ... + + pi.*(Alpha-1) .* cos(pi.*(Alpha-1)./(4.*Alpha)) ... + ); +end + +ind = 1:length(n); +ind([ind1 ind2]) = []; +nind = n(ind); + +b(ind) = -4.*Alpha./fa .* ( cos((1+Alpha).*2.*pi.*nind./Tsym) + ... + sin((1-Alpha).*2.*pi.*nind./Tsym) ./ (8.*Alpha.*nind./Tsym) ) ... + ./ (pi .* sqrt(1./(2./Tsym)) .* ((8.*Alpha.*nind./Tsym).^2 - 1)); + +b = sqrt(2./Tsym) .* b; +k=0.5*fa*Tsym; diff --git a/common/radio/fir_rrc_bak.m b/common/radio/fir_rrc_bak.m new file mode 100755 index 0000000..b1f0c23 --- /dev/null +++ b/common/radio/fir_rrc_bak.m @@ -0,0 +1,26 @@ +% function b = fir_rrc(fa, T, N, alpha) +% Root Raised Cosine FIR-Filter + +function b = fir_rrc(fa, T, N, alpha) + +n = -(N-1)/2:1:(N-1)/2; +t = n/fa; + +A = 4*alpha/(pi*sqrt(T)) + +nom = 4*alpha*cos((1+alpha)*pi*t/T) + T*sin((1-alpha)*pi*t/T); +denom = pi*sqrt(T)*(1 - (4*alpha*t/T).^2); + +% Find zero elements +ze = find(abs(alpha*t/T) == 0.25) + +%nom(ze) = sqrt(T)*(1 - alpha*(1 - (4/pi))); +%denom(ze) = 1; + +b = nom./denom; + + + + + + diff --git a/common/rgb2xyz.m b/common/rgb2xyz.m new file mode 100755 index 0000000..e22dc40 --- /dev/null +++ b/common/rgb2xyz.m @@ -0,0 +1,44 @@ +% function [xyz] = rgb2xyz(rgb) +% +% // C-Code +% if ( var_R > 0.04045 ) var_R = ( ( var_R + 0.055 ) / 1.055 ) ^ 2.4 +% else var_R = var_R / 12.92 +% if ( var_G > 0.04045 ) var_G = ( ( var_G + 0.055 ) / 1.055 ) ^ 2.4 +% else var_G = var_G / 12.92 +% if ( var_B > 0.04045 ) var_B = ( ( var_B + 0.055 ) / 1.055 ) ^ 2.4 +% else var_B = var_B / 12.92 +% +% var_R = var_R * 100 +% var_G = var_G * 100 +% var_B = var_B * 100 +% +% // Observer. = 2°, Illuminant = D65 +% X = var_R * 0.4124 + var_G * 0.3576 + var_B * 0.1805 +% Y = var_R * 0.2126 + var_G * 0.7152 + var_B * 0.0722 +% Z = var_R * 0.0193 + var_G * 0.1192 + var_B * 0.9505 +% +% Source http://www.easyrgb.com + +function [xyz] = rgb2xyz(rgb) +[m,n] = size(rgb); + +if ((m ~= 3) & (n == 3)) + rgb = rgb'; +end; + +% Observer = 2°, Illuminant = D65 +t_mat = [[ 0.412453 0.357580 0.180423 ] + [ 0.212671 0.715160 0.072169 ] + [ 0.019334 0.119193 0.950227 ]]; + +rgb_k = f1(rgb); +xyz = t_mat * rgb_k; + +function rgb_k = f1(rgb) +thresh_n = 0.04045 ; + +[idx_m, idx_n] = find(rgb > thresh_n); +rgb_k(idx_m,idx_n) = ((rgb(idx_m,idx_n) + 0.055)/1.055).^(2.4); + +[idx_m, idx_n] = find(rgb <= thresh_n); +rgb_k(idx_m,idx_n) = rgb(idx_m,idx_n)/12.92; diff --git a/common/rls.m b/common/rls.m new file mode 100755 index 0000000..9cf729e --- /dev/null +++ b/common/rls.m @@ -0,0 +1,55 @@ +% RLS-Algorithmus +% --------------- +% function [dd,e,wf,fpo] = rls(x,d,N,mu,rho,wi) +% +% Parameter: +% x : Eingangssignal +% d : erwünschtes Signal (Referenzsignal) +% N : Filterordnung +% mu : konstante Schrittweite +% rho : Vergessensfaktor +% wi : Startwerte der Filterkoeffizienten im Zeitbereich +% +% Rückgabewerte: +% dd : Schätzung d' des Referenzsignals d aus x +% e : Fehler d - d' +% wf : Filterkoeffizienten am Ende der Adaption im Zeitbereich +% fpo : Anzahl der benötigten Floating-Point Operationen + +% ------------------------------------------------------------------------- +% Datum : 27.3.2002 +% Autor : Jens Ahrensfeld +% Thema : Diplomarbeit +% Datei : rls.m +% +% ------------------------------------------------------------------------- +function [dd,e,wf,fpo] = rls(x,d,N,mu,rho,wi) + +% ------------------------------------------------------------------------- +% Initialisierung RLS +% ------------------------------------------------------------------------- +wf = wi; +dd = zeros(length(x),1); +e = zeros(length(x),1); +eta = 1000000; +R_ = eta * eye(N); % inverse Autokorrelationsmatrix R +flops(0); + +% ------------------------------------------------------------------------- +% RLS Adaptation +% ------------------------------------------------------------------------- +for n = N : length(x) + xs = x(n:-1:n-N+1); % Eingangsvektor + dd(n) = xs' * wf; % Filterung + e(n) = d(n) - dd(n); % a priori-Fehler + z = R_ * xs; % gefilterter Datenvektor z + vn = 1/(rho + xs'*z); % Normierungskonstante + zn = vn*z; % Normierung + wf = wf + mu*e(n)*zn; % Aktualisierung der Filterkoeffizienten + R_ = 1/rho * (R_ - zn*xs'*R_); % Aktualisierung der inversen + % Autokorrelationsmatrix R +end ; +fpo = flops; + +% ------------------------------------------------------------------------- +% Ende rls.m diff --git a/common/rls1.m b/common/rls1.m new file mode 100755 index 0000000..f1b1ad1 --- /dev/null +++ b/common/rls1.m @@ -0,0 +1,60 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% RLS Algorithm % +% % +% Written By: Sundar Sankaran and A. A. (Louis) Beex % +% DSP Research Laboratory % +% Dept. of Electrical and Comp. Engg % +% Virginia Tech % +% Blacksburg VA 24061-0111 % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +randn('seed', 0) ; +rand('seed', 0) ; + +NoOfData = 8000 ; % Set no of data points used for training +Order = 32 ; % Set the adaptive filter order + +Lambda = 0.98 ; % Set the forgetting factor +Delta = 0.001 ; % R initialized to Delta*I + +x = randn(NoOfData, 1) ;% Input assumed to be white +h = rand(Order, 1) ; % System picked randomly +d = filter(h, 1, x) ; % Generate output (desired signal) + +% Initialize RLS + +P = Delta * eye ( Order, Order ) ; +w = zeros ( Order, 1 ) ; + +% RLS Adaptation + +for n = Order : NoOfData ; + + u = x(n:-1:n-Order+1) ; + pi_ = u' * P ; + k = Lambda + pi_ * u ; + K = pi_'/k; + e(n) = d(n) - w' * u ; + w = w + K * e(n) ; + PPrime = K * pi_ ; + P = ( P - PPrime ) / Lambda ; + w_err(n) = norm(h - w) ; + +end ; + + +% Plot results + +figure ; +plot(20*log10(abs(e))) ; +title('Learning Curve') ; +xlabel('Iteration Number') ; +ylabel('Output Estimation Error in dB') ; + +figure ; +semilogy(w_err) ; +title('Weight Estimation Error') ; +xlabel('Iteration Number') ; +ylabel('Weight Error in dB') ; diff --git a/common/rls2.m b/common/rls2.m new file mode 100755 index 0000000..48907dc --- /dev/null +++ b/common/rls2.m @@ -0,0 +1,71 @@ +%=================================================================== +%------------------------------------------------------------------- +% +% Adaptive Filter +% +% Simulationen +% +% [E,W,w,inv_R]=rls(N,X,D,w_start,rho) +% +% RLS-Algorithmus +% +% E = Fehlersignal E[.] +% W = Filterkoeffizienten im zeitlichen Verlauf +% w = Filterkoeffizienten am Ende der Adaption +% inv_R = Inverse deterministische Korrelationsmatrix +% nach Adaptionsende +% N = Anz. Filterkoeffizienten = Filterordnung+1 +% X = Filtereingang X[.] +% D = erwuenschtes Signal D[.] +% w_start = Startwerte Filterkoeffizienten +% rho = Vergessensfaktor +% +%------------------------------------------------------------------- + +% +% author: Markus Hofbauer +% ISI, ETH Zuerich (Switzerland) +% +% created: 7/2000 +% +% +%------------------------------------------------------------------- +% +% File : rls.m +% +% Startfile: simX.m +% +%------------------------------------------------------------------- +%=================================================================== + +function [E,W,w,inv_R]=rls(N,X,D,w_start,rho) + +p0=1000000; % Initialisierung von +inv_R=p0*eye(N); % inv_R + +adaptlen=length(X); +w= w_start; +W=zeros(N,adaptlen); +E=zeros(adaptlen,1); + +%------------------------------------------------------------------- +% RLS Update Loop +%------------------------------------------------------------------- + +for i=N:adaptlen + + W(:,i)=w; + + x = X(i:-1:i-N+1); % Eingangsvektor (x[k],x[k-1],..,x[k-N+1]) + y=x'*w; % Filterausgang + e=D(i)-y; % Fehler + + c=1/(rho+x'*inv_R*x); + inv_R=1/rho*(inv_R-c*inv_R*x*x'*inv_R); % Aufdatierung von inv_R + w=w+inv_R*e*x; % Aufdatierung von w + + E(i)=e; +end; + + + \ No newline at end of file diff --git a/common/rms.m b/common/rms.m new file mode 100755 index 0000000..09d193f --- /dev/null +++ b/common/rms.m @@ -0,0 +1,13 @@ +% Power estimation on discrete signals +function y = rms(x,Nw) + +N = length(x); +y = zeros(1,N); + +y(1) = x(1)^2; + +for k = 2:N-Nw, + for n = k:Nw+k, + y(n) = y(n-1)*(1-1/Nw) + 1/Nw*x(n)^2; + end; +end; \ No newline at end of file diff --git a/common/sidex.m b/common/sidex.m new file mode 100755 index 0000000..204fd41 --- /dev/null +++ b/common/sidex.m @@ -0,0 +1,33 @@ +% sidex.m - Demonstration of using FFT cross-correlation to compute +% the impulse response of a filter given its input and output. +% This is called "FIR system identification". + +Nx = 256; % input signal length +Nh = 100; % filter length +Ny = Nx+Nh-1; % max output signal length + +% FFT size to accommodate cross-correlation: +Nfft = 2^nextpow2(Nx+Ny-1); % want power of 2 for FFT + +%x = rand(1,Nx); % input signal = noise +x = 1:Nx; % input signal = ramp +h = [1:Nh]; % the filter +xzp = [x,zeros(1,Nfft-Nx)]; % zero-padded input signal +yzp = filter(h,1,xzp); % apply the filter +X = fft(xzp); % input spectrum +Y = fft(yzp); % output spectrum +Rxx = conj(X) .* X; % energy spectrum of x +Rxy = conj(X) .* Y; % cross-energy spectrum of x and y +Hxy = Rxy ./ Rxx; % should be the freq. response +hxy = ifft(Hxy); % should be the imp. response + +hxy(1:Nh) % print estimated impulse response +%freqz(hxy,1,Nfft); % plot estimated frequency response +plot(1:lge(hxy),real(hxy)); + +err = norm(hxy - [h,zeros(1,Nfft-Nh)])/norm(h); +disp(sprintf('Impulse Response Error = %0.14f%%',100*err)); + +err = norm(Hxy - fft([h,zeros(1,Nfft-Nh)]))/norm(h); +disp(sprintf('Frequency Response Error = %0.14f%%',100*err)); + diff --git a/common/sm.m b/common/sm.m new file mode 100755 index 0000000..6896768 --- /dev/null +++ b/common/sm.m @@ -0,0 +1,21 @@ +% Power estimation on discrete signals +function [y,yf] = rms(x,Nw,yi) + +N = length(x); +y = zeros(1,N); + + +Nb = ceil(N/Nw); +xp = zeros(1,Nb*Nw); +yp = zeros(1,Nb*Nw); + +xp(1:N) = x(1:N); +yp(1) = abs(x(1)); + +for k = 1:Nb-1, + for n = (k-1)*Nw+2:k*Nw+2, + yp(n) = yp(n-1)*(1-1/Nw) + 1/Nw*abs(x(n)); + end; +end; + +y = yp(1:N); \ No newline at end of file diff --git a/common/sm2.m b/common/sm2.m new file mode 100755 index 0000000..b8a0bfd --- /dev/null +++ b/common/sm2.m @@ -0,0 +1,9 @@ +% Power estimation on discrete signals +function [y,yf] = sm2(x,yi,arf) + +N = length(x); +x(1) = x(1) + yi*(1/arf-1); +x = abs(x); +a = [1,-(1-arf)]; +y = filter(arf,a,x); +yf = y(N); diff --git a/common/smhs.m b/common/smhs.m new file mode 100755 index 0000000..ade6c34 --- /dev/null +++ b/common/smhs.m @@ -0,0 +1,21 @@ +% Power estimation on discrete signals +function [y,yf] = smhs(x,Nw1,Nw2,yi) + +N = length(x); +y = zeros(1,N); +x = abs(x); + +if (x(1) > yi) + y(1) = (1-1/Nw1)*yi + (1/Nw1)*x(1); +else + y(1) = (1-1/Nw2)*yi + (1/Nw2)*x(1); +end; + +for n = 2:N, + if (x(n) > y(n-1)) + y(n) = (1-1/Nw1)*y(n-1) + (1/Nw1)*x(n); + else + y(n) = (1-1/Nw2)*y(n-1) + (1/Nw2)*x(n); + end; +end; +yf = y(N); diff --git a/common/testcoef.m b/common/testcoef.m new file mode 100755 index 0000000..19868e6 --- /dev/null +++ b/common/testcoef.m @@ -0,0 +1,10 @@ +% testcoef(N,k) +function b = testcoef(N,K) + +b = ones(N,1); + +k = 1; +for m=1:N, + b(m) = b(m)*k*(-1)^(m+1); + k = K*k; +end; diff --git a/common/tristimulus.m b/common/tristimulus.m new file mode 100755 index 0000000..50a85a3 --- /dev/null +++ b/common/tristimulus.m @@ -0,0 +1,73 @@ +% function [ref] = tristimulus(illuminant, observer) +% + +function [ref] = tristimulus(illuminant, observer) + +ref = [0.95047 1 1.08883]'; + +switch observer + case {'2°'} + switch illuminant + case {'A'} + ref = [1.09850 1 0.35585]'; + case {'C'} + ref = [0.98074 1 1.18232]'; + case {'D50'} + ref = [0.96422 1 0.82521]'; + case {'D55'} + ref = [0.95682 1 0.92149]'; + case {'D65'} + ref = [0.95047 1 1.08883]'; + case {'D75'} + ref = [0.94972 1 1.22638]'; + case {'F2'} + ref = [0.99187 1 0.67395]'; + case {'F7'} + ref = [0.95044 1 1.08755]'; + case {'F11'} + ref = [1.00966 1 0.64370]'; + + case {'Incandescent'} % same as 'A' + ref = [1.09850 1 0.35585]'; + case {'Daylight'} % same as 'D65' + ref = [0.95047 1 1.08883]'; + case {'Fluorescent'} % same as 'F2' + ref = [0.99187 1 0.67395]'; + + otherwise + ref = [0.95047 1 1.08883]'; + end; + + case {'10°'} + switch illuminant + case {'A'} + ref = [1.11144 1 0.35200]'; + case {'C'} + ref = [0.97285 1 1.16145]'; + case {'D50'} + ref = [0.96720 1 0.81427]'; + case {'D55'} + ref = [0.95799 1 0.90926]'; + case {'D65'} + ref = [0.94811 1 1.07304]'; + case {'D75'} + ref = [0.94416 1 1.20641]'; + case {'F2'} + ref = [1.03280 1 0.69026]'; + case {'F7'} + ref = [0.95792 1 1.07687]'; + case {'F11'} + ref = [1.03866 1 0.65627]'; + + case {'Incandescent'} % same as 'A' + ref = [1.11144 1 0.35200]'; + case {'Daylight'} % same as 'D65' + ref = [0.94811 1 1.07304]'; + case {'Fluorescent'} % same as 'F2' + ref = [1.03280 1 0.69026]'; + + otherwise + ref = [0.94811 1 1.07304]'; + end; +end; + diff --git a/common/undcc.m b/common/undcc.m new file mode 100755 index 0000000..5416448 --- /dev/null +++ b/common/undcc.m @@ -0,0 +1,27 @@ +% Erzeugt Filterkoeffizienten für FIR-Filter zur Eliminierung des Gleichanteils +% ----------------------------------------------------------------------------- +% +% Aufruf +% Function [coeff] = undcc(N) +% +% Parameter: +% N : Bestimmt Grösse des Koeffizientenvektors +% +% Rückgabewerte: +% coeff : Koeffizientenvektor, Länge N + +% ----------------------------------------------------------------------------- +% Datum : 26.07.2001 +% Autor : Jens Ahrensfeld +% Thema : Diplomarbeit +% Datei : undcc.m +% Benötigte Dateien: +% +% ----------------------------------------------------------------------------- +function [coeff] = undcc(N) + +% Ansatz: y = x - arith. Mittelwert(x) +coeff = [1.0-1/N ; -1/N*ones(N-1,1)]; + +% ----------------------------------------------------------------------------- +% Ende undcc.m \ No newline at end of file diff --git a/common/wav2dat.m b/common/wav2dat.m new file mode 100755 index 0000000..1582bee --- /dev/null +++ b/common/wav2dat.m @@ -0,0 +1,16 @@ +% function y = wav2dat(name) +% + +function y = wav2dat(name) + +datfile = sprintf('%s.dat', name); +wavfile = sprintf('%s.wav', name); + +[data, fa, fmt] = wavr16(wavfile); +nBits = fmt(6); +data_float = 2^-(nBits-1)*data; + +fid = fopen(datfile, 'w'); +fwrite(fid, data_float, 'float32'); +fclose(fid); + diff --git a/common/wavr16.m b/common/wavr16.m new file mode 100755 index 0000000..0f5354c --- /dev/null +++ b/common/wavr16.m @@ -0,0 +1,67 @@ +function [y,Fs,Format]=wavr16(wavefile) +%WAVREAD Load Microsoft Windows 3.1 .WAV format sound files. +% [y]=WAVREAD(wavefile) loads a .WAV format file specified by "wavefile", +% returning the sampled data in variable "y". The .WAV extension +% in the filename is optional. +% +% [y,Fs]=WAVREAD(wavefile) loads a .WAV format file specified by +% "wavefile", returning the sampled data in variable "y" and the +% sample rate in variable "Fs". +% +% [y,Fs,Format]=WAVREAD(wavefile) loads a .WAV format file specified by +% "wavefile",returning the sampled data in variable "y", the sample +% rate in variable "Fs", and the .WAV file format information in +% variable "Format". The format information is returned as a 6 element +% vector with the following order: +% +% Format(1) Data format (always PCM) +% Format(2) Number of channels +% Format(3) Sample Rate (Fs) +% Format(4) Average bytes per second (sampled) +% Format(5) Block alignment of data +% Format(6) Bits per sample +% +% Note: WAVREAD currently supports only 8-bit single channel data. +% +% See also WAVWRITE. + +% Copyright (c) 1984-94 by The MathWorks, Inc. + +if nargin~=1 + error('WAVREAD takes one argument, which is the name of the .WAV file'); +end + +if isempty(findstr(wavefile,'.')) + wavefile=[wavefile,'.wav']; +end + +fid=fopen(wavefile,'rb','l'); +if fid ~= -1 + % read riff chunk + header=fread(fid,4,'uchar'); + header=fread(fid,1,'ulong'); + header=fread(fid,4,'uchar'); + + % read format sub-chunk + header=fread(fid,4,'uchar'); + header=fread(fid,1,'ulong'); + + Format(1)=fread(fid,1,'ushort'); % PCM format + Format(2)=fread(fid,1,'ushort'); % 1 channel + Fs=fread(fid,1,'ulong'); % samples per second + Format(3)=Fs; + Format(4)=fread(fid,1,'ulong'); % average bytes per second + Format(5)=fread(fid,1,'ushort'); % block alignment + Format(6)=fread(fid,1,'ushort'); % bits per sample + + + % read data sub-chunck + header=fread(fid,4,'uchar'); + nsamples=fread(fid,1,'ulong'); + y=fread(fid,nsamples,'short'); + fclose(fid); +end + +if fid == -1 + error('Can''t open .WAV file for input!'); +end; diff --git a/common/wavw16.m b/common/wavw16.m new file mode 100755 index 0000000..7ae8875 --- /dev/null +++ b/common/wavw16.m @@ -0,0 +1,70 @@ +function wavw16(waveData,sRate,wavefile) +%WAVWRITE Saves Microsoft Windows 3.1 .WAV format sound files. +% WAVWRITE(y,Fs,wavefile) saves a .WAV format file specified by "wavefile". +% +% The input arguments for WAVWRITE are as follows: +% +% y The sampled data to save (8 bit max) +% Fs The rate at which the data was sampled +% wavefile A string containing the name of the .WAV file to create +% +% Note: WAVWRITE will create an 8-bit, simgle channel wave file. Non 8-bit +% sample data will be truncated. +% +% See also WAVREAD. + +% Copyright (c) 1984-94 by The MathWorks, Inc. + +if nargin~=3 + error('WAVWRITE needs three arguments!'); +end +if isstr(waveData) %old symtax, reorder args + tmp = waveData; + waveData = wavefile; + wavefile = sRate; + sRate = tmp; +end + +if isempty(findstr(wavefile,'.')) + wavefile=[wavefile,'.wav']; +end + +fid=fopen(wavefile,'wb','l'); + +if fid ~= -1 + [m,n]=size(waveData); + nsamples=m*n; + + riffsize=36+nsamples*2; + + % write riff chunk + fwrite(fid,'RIFF','uchar'); + fwrite(fid,riffsize,'ulong'); + fwrite(fid,'WAVE','uchar'); + + % write format sub-chunk + fwrite(fid,'fmt ','uchar'); + fwrite(fid,16,'ulong'); + + fwrite(fid,1,'ushort'); % PCM format + fwrite(fid,1,'ushort'); % 1 channel + fwrite(fid,sRate,'ulong'); % samples per second + fwrite(fid,2*sRate,'ulong'); % average bytes per second + fwrite(fid,2,'ushort'); % block alignment + fwrite(fid,16,'ushort'); % bits per sample + + + % write data sub-chunck + fwrite(fid,'data','uchar'); + fwrite(fid,2*nsamples,'ulong'); + fwrite(fid,waveData,'short'); + + if ((ceil(nsamples/2))~=(nsamples/2)) + fwrite(fid,0,'short'); + end; + fclose(fid); +end; + +if fid == -1 + error('Can''t open .WAV file for output!'); +end; diff --git a/common/whamm.m b/common/whamm.m new file mode 100755 index 0000000..ff7a609 --- /dev/null +++ b/common/whamm.m @@ -0,0 +1,22 @@ +% Hamming Window +% ------------------------------------ +% +% Aufruf: +% Function [y]=whamm(x) +% +% Parameter: +% +% Rückgabewerte: + +% ------------------------------------------------------------------------- +% Datum : ..2001 +% Autor : Jens Ahrensfeld +% Thema : Diplomarbeit +% Datei : newfunc.m +% Benötigte Dateien: +% +% ------------------------------------------------------------------------- +function y=whamm(x) + +n = 0:lge(x)-1; +y = x.*(0.54-0.46*cos(2*pi*n/(lge(x)-1)))'; diff --git a/common/winhamm.m b/common/winhamm.m new file mode 100755 index 0000000..a18829e --- /dev/null +++ b/common/winhamm.m @@ -0,0 +1,22 @@ +% Hamming Window +% ------------------------------------ +% +% Aufruf: +% Function [y]=whamm(x) +% +% Parameter: +% +% Rückgabewerte: + +% ------------------------------------------------------------------------- +% Datum : ..2001 +% Autor : Jens Ahrensfeld +% Thema : Diplomarbeit +% Datei : newfunc.m +% Benötigte Dateien: +% +% ------------------------------------------------------------------------- +function y=winhamm(x) + +n = 0:lge(x)-1; +y = x.*(0.54-0.46*cos(2*pi*n/(lge(x)-1))); diff --git a/common/wr16.m b/common/wr16.m new file mode 100755 index 0000000..eb505e1 --- /dev/null +++ b/common/wr16.m @@ -0,0 +1,67 @@ +function [y,Fs,Format]=wr16(wavefile) +%WAVREAD Load Microsoft Windows 3.1 .WAV format sound files. +% [y]=WAVREAD(wavefile) loads a .WAV format file specified by "wavefile", +% returning the sampled data in variable "y". The .WAV extension +% in the filename is optional. +% +% [y,Fs]=WAVREAD(wavefile) loads a .WAV format file specified by +% "wavefile", returning the sampled data in variable "y" and the +% sample rate in variable "Fs". +% +% [y,Fs,Format]=WAVREAD(wavefile) loads a .WAV format file specified by +% "wavefile",returning the sampled data in variable "y", the sample +% rate in variable "Fs", and the .WAV file format information in +% variable "Format". The format information is returned as a 6 element +% vector with the following order: +% +% Format(1) Data format (always PCM) +% Format(2) Number of channels +% Format(3) Sample Rate (Fs) +% Format(4) Average bytes per second (sampled) +% Format(5) Block alignment of data +% Format(6) Bits per sample +% +% Note: WAVREAD currently supports only 8-bit single channel data. +% +% See also WAVWRITE. + +% Copyright (c) 1984-94 by The MathWorks, Inc. + +if nargin~=1 + error('WAVREAD takes one argument, which is the name of the .WAV file'); +end + +if findstr(wavefile,'.')==[ ] + wavefile=[wavefile,'.wav']; +end + +fid=fopen(wavefile,'rb','l'); +if fid ~= -1 + % read riff chunk + header=fread(fid,4,'uchar'); + header=fread(fid,1,'ulong'); + header=fread(fid,4,'uchar'); + + % read format sub-chunk + header=fread(fid,4,'uchar'); + header=fread(fid,1,'ulong'); + + Format(1)=fread(fid,1,'ushort'); % PCM format + Format(2)=fread(fid,1,'ushort'); % 1 channel + Fs=fread(fid,1,'ulong'); % samples per second + Format(3)=Fs; + Format(4)=fread(fid,1,'ulong'); % average bytes per second + Format(5)=fread(fid,1,'ushort'); % block alignment + Format(6)=fread(fid,1,'ushort'); % bits per sample + + + % read data sub-chunck + header=fread(fid,4,'uchar'); + nsamples=fread(fid,1,'ulong'); + y=fread(fid,nsamples,'short'); + fclose(fid); +end + +if fid == -1 + error('Can''t open .WAV file for input!'); +end; diff --git a/common/xlte.m b/common/xlte.m new file mode 100755 index 0000000..c217800 --- /dev/null +++ b/common/xlte.m @@ -0,0 +1,41 @@ +% Long-Term Magnitude Estimation von x +% ------------------------------------- +% +% Aufruf: +% Function [xl,xlf] = xlte(x,xs,af,pl,xli) +% +% Parameter: +% x : +% ar : +% af : +% xsi : +% +% Rückgabewerte: +% xs : +% xsf : + +% ------------------------------------------------------------------------- +% Datum : 26.07.2001 +% Autor : Jens Ahrensfeld +% Thema : Diplomarbeit +% Datei : xste.m +% Benötigte Dateien: +% +% ------------------------------------------------------------------------- +function [xl,xlf] = xlte(x,xs,af,pl,xli) + +N = lge(x); +x = abs(x); +xl(1) = xli; + +for n = 2:N, + if(xs(n) > pl*xl(n-1)) + xl(n) = xl(n-1); + else + xl(n) = af*x(n) + (1-af)*xl(n-1); + end; +end; +xlf = xl(N); + +% ------------------------------------------------------------------------- +% Ende xlte.m diff --git a/common/xste.m b/common/xste.m new file mode 100755 index 0000000..810d094 --- /dev/null +++ b/common/xste.m @@ -0,0 +1,41 @@ +% Short-Term Magnitude Estimation von x +% ------------------------------------- +% +% Aufruf: +% Function [xs,xsf] = xste(x,ar,af,xsi) +% +% Parameter: +% x : +% ar : +% af : +% xsi : +% +% Rückgabewerte: +% xs : +% xsf : + +% ------------------------------------------------------------------------- +% Datum : 26.07.2001 +% Autor : Jens Ahrensfeld +% Thema : Diplomarbeit +% Datei : xste.m +% Benötigte Dateien: +% +% ------------------------------------------------------------------------- +function [xs,xsf] = xste(x,ar,af,xsi) + +N = lge(x); +x = abs(x); +xs(1) = xsi; + +for n = 2:N, + if(x(n) > xs(n-1)) + xs(n) = ar*x(n) + (1-ar)*xs(n-1); + else + xs(n) = af*x(n) + (1-af)*xs(n-1); + end; +end; +xsf = xs(N); + +% ------------------------------------------------------------------------- +% Ende xste.m diff --git a/common/xyz2cielab.m b/common/xyz2cielab.m new file mode 100755 index 0000000..10deefe --- /dev/null +++ b/common/xyz2cielab.m @@ -0,0 +1,47 @@ +% function [lab] = xyz2cielab(xyz, ref_xyz) +% +% // C-Code +% var_X = X / 95.047 //Observer = 2°, Illuminant = D65 +% var_Y = Y / 100.000 +% var_Z = Z / 108.883 +% +% if ( var_X > 0.008856 ) var_X = var_X ^ ( 1/3 ) +% else var_X = ( 7.787 * var_X ) + ( 16 / 116 ) +% if ( var_Y > 0.008856 ) var_Y = var_Y ^ ( 1/3 ) +% else var_Y = ( 7.787 * var_Y ) + ( 16 / 116 ) +% if ( var_Z > 0.008856 ) var_Z = var_Z ^ ( 1/3 ) +% else var_Z = ( 7.787 * var_Z ) + ( 16 / 116 ) +% +% CIE-L* = ( 116 * var_Y ) - 16 +% CIE-a* = 500 * ( var_X - var_Y ) +% CIE-b* = 200 * ( var_Y - var_Z ) +% +% Source http://www.easyrgb.com + +function [lab] = xyz2cielab(xyz, ref_xyz) + +[m,n] = size(xyz); +if ((m ~= 3) & (n == 3)) + xyz = xyz' +end; + +xyz(1,:) = xyz(1,:)./ref_xyz(1,:); +xyz(2,:) = xyz(2,:)./ref_xyz(2,:); +xyz(3,:) = xyz(3,:)./ref_xyz(3,:); +xyz_k = f1(xyz); + +% L +lab(1,:) = 116*xyz_k(2,:) - 16; +% a +lab(2,:) = 500*(xyz_k(1,:) - xyz_k(2,:)); +% b +lab(3,:) = 200*(xyz_k(2,:) - xyz_k(3,:)); + +function xyz_k = f1(xyz) +thresh_n = 0.008856; + +[idx_m, idx_n] = find(xyz > thresh_n); +xyz_k(idx_m,idx_n) = xyz(idx_m,idx_n).^(1/3); + +[idx_m, idx_n] = find(xyz <= thresh_n); +xyz_k(idx_m,idx_n) = 7.787 * xyz(idx_m,idx_n) + 16/116; diff --git a/common/xyz2rgb.m b/common/xyz2rgb.m new file mode 100755 index 0000000..7c3de33 --- /dev/null +++ b/common/xyz2rgb.m @@ -0,0 +1,47 @@ +% function [rgb] = xyz2rgb(xyz) +% +% // C-Code +% ref_X = 95.047 //Observer = 2°, Illuminant = D65 +% ref_Y = 100.000 +% ref_Z = 108.883 +% +% var_X = X / 100 //X = From 0 to ref_X +% var_Y = Y / 100 //Y = From 0 to ref_Y +% var_Z = Z / 100 //Z = From 0 to ref_Y +% +% var_R = var_X * 3.2406 + var_Y * -1.5372 + var_Z * -0.4986 +% var_G = var_X * -0.9689 + var_Y * 1.8758 + var_Z * 0.0415 +% var_B = var_X * 0.0557 + var_Y * -0.2040 + var_Z * 1.0570 +% +% if ( var_R > 0.0031308 ) var_R = 1.055 * ( var_R ^ ( 1 / 2.4 ) ) - 0.055 +% else var_R = 12.92 * var_R +% if ( var_G > 0.0031308 ) var_G = 1.055 * ( var_G ^ ( 1 / 2.4 ) ) - 0.055 +% else var_G = 12.92 * var_G +% if ( var_B > 0.0031308 ) var_B = 1.055 * ( var_B ^ ( 1 / 2.4 ) ) - 0.055 +% else var_B = 12.92 * var_B +% +% Source http://www.easyrgb.com + +function [rgb] = xyz2rgb(xyz) +[m,n] = size(xyz); + +if ((m ~= 3) & (n == 3)) + xyz = xyz'; +end; + +% Observer = 2°, Illuminant = D65 +t_mat = [[ 3.240479 -1.537150 -0.498535 ] + [ -0.969256 1.875992 0.041556 ] + [ 0.055648 -0.204043 1.057311 ]]; + +rgb_k = t_mat * xyz; +rgb = f1(rgb_k); + +function rgb = f1(rgb_k) +thresh_n = 0.0031308; + +[idx_m,idx_n] = find(rgb_k > thresh_n); +rgb(idx_m,idx_n) = 1.055*(rgb_k(idx_m,idx_n).^(1/2.4)) - 0.055; + +[idx_m,idx_n] = find(rgb_k <= thresh_n); +rgb(idx_m,idx_n) = 12.92*rgb_k(idx_m,idx_n);