From 1f5b0b0570e142ad12c58e50f0e8e9b642951719 Mon Sep 17 00:00:00 2001 From: Jens Ahrensfeld Date: Thu, 23 Jul 2020 11:10:12 +0000 Subject: [PATCH] - added git-svn-id: http://moon:8086/svn/matlab/trunk@128 801c6759-fa7c-4059-a304-17956f83a07c --- equalizer/boost.m | 28 ++++++++++++++++ equalizer/boost2.m | 29 ++++++++++++++++ equalizer/boost3.m | 48 +++++++++++++++++++++++++++ equalizer/boost3a.m | 43 ++++++++++++++++++++++++ equalizer/eq_eval.m | 81 +++++++++++++++++++++++++++++++++++++++++++++ kalman_eval.m | 2 +- stock/diba_aktie.m | 10 ++++++ tempsensor_eval.m | 48 +++++++++++++++++++++++++++ 8 files changed, 288 insertions(+), 1 deletion(-) create mode 100644 equalizer/boost.m create mode 100644 equalizer/boost2.m create mode 100644 equalizer/boost3.m create mode 100644 equalizer/boost3a.m create mode 100644 equalizer/eq_eval.m create mode 100644 stock/diba_aktie.m create mode 100644 tempsensor_eval.m diff --git a/equalizer/boost.m b/equalizer/boost.m new file mode 100644 index 0000000..d870d19 --- /dev/null +++ b/equalizer/boost.m @@ -0,0 +1,28 @@ +function [B,A] = boost(fs, fc, bw, gain, P); +%BOOST - Design a boost filter at given gain, center +% frequency fc, bandwidth bw, and sampling rate fs +% (default = 1). +% +% J.O. Smith 11/28/02 +% Reference: Zolzer: Digital Audio Signal Processing, p. 124 + +Q = fs/bw; +wcT = 2*pi*fc/fs; + +K=tan(wcT/2); +V=gain; + +b0 = 1 + V*K/Q + K^2; +b1 = 2*(K^2 - 1); +b2 = 1 - V*K/Q + K^2; +a0 = 1 + K/Q + K^2; +a1 = 2*(K^2 - 1); +a2 = 1 - K/Q + K^2; +A(1,:) = [a0 a1 a2] /a0; +B(1,:) = [b0 b1 b2] /a0; + +if nargout==0 + freqz(B(1,:),A(1,:)); + title('Boost Frequency Response') +end + diff --git a/equalizer/boost2.m b/equalizer/boost2.m new file mode 100644 index 0000000..ba7e512 --- /dev/null +++ b/equalizer/boost2.m @@ -0,0 +1,29 @@ +function [B,A] = boost2(fs, fc, bw, gain); +%BOOST - Design a boost filter at given gain, center +% frequency fc, bandwidth bw, and sampling rate fs +% (default = 1). +% +% Cookbook formulae for audio EQ biquad filter coefficients +% by Robert Bristow-Johnson + +Q = fs/bw; +w0 = 2*pi*fc/fs; + +A=sqrt(gain); +alpha = sin(w0)/(2*Q) + +b0 = 1 + alpha*A +b1 = -2*cos(w0) +b2 = 1 - alpha*A +a0 = 1 + alpha/A +a1 = -2*cos(w0) +a2 = 1 - alpha/A + +A = [a0 a1 a2] / a0; +B = [b0 b1 b2] / a0; + +if nargout==0 + freqz(B,A); + title('Boost Frequency Response') +end + diff --git a/equalizer/boost3.m b/equalizer/boost3.m new file mode 100644 index 0000000..5b92b47 --- /dev/null +++ b/equalizer/boost3.m @@ -0,0 +1,48 @@ +function [Bp,Ap] = boost3(fs, fc, bw, gain, P); +%BOOST - Design a boost filter at given gain, center +% frequency fc, bandwidth bw, and sampling rate fs +% (default = 1). +% +% Cookbook formulae for audio EQ biquad filter coefficients +% by Robert Bristow-Johnson + +Q = bw/fs; + +Bp = zeros(P, 3); +Ap = zeros(P, 3); +for p=1:P, + kp = IIRCalcQp(p, P); + [Bp(p,:),Ap(p,:)] = peaking(fs, fc, kp*Q, gain); +endfor + +B = Bp(1,:); +A = Ap(1,:); +if nargout==0 + close + freqz(B,A); + title('Boost Frequency Response') +end + +endfunction + +function [B,A] = peaking(fs, fc, Q, K) + k1=-cos(2*pi*fc/fs); + + if K < 1 + k2=(1-tan(Q/K/2))/(1+tan(Q/K/2)); + else + k2=(1-tan(Q/2))/(1+tan(Q/2)); + end + Pz = [1 k1*(1+k2) k2]; % define denominator coefficients + Qz = [k2 k1*(1+k2) 1]; % define numerator coefficients + Num = (Pz*(1+K) + Qz*(1-K))/2; + Den = Pz; + B = Num; + A = Den; + +endfunction + +function kp = IIRCalcQp(p, P) + kp = 1.0/(2*sin(pi*(2*p-1)/(2*P))); + +endfunction diff --git a/equalizer/boost3a.m b/equalizer/boost3a.m new file mode 100644 index 0000000..db393bf --- /dev/null +++ b/equalizer/boost3a.m @@ -0,0 +1,43 @@ +function [Bp,Ap] = boost3a(fs, fc, bw, gain, P); +%BOOST - Design a boost filter at given gain, center +% frequency fc, bandwidth bw, and sampling rate fs +% (default = 1). +% +% Cookbook formulae for audio EQ biquad filter coefficients +% by Robert Bristow-Johnson + +Q = bw/fs; + +Bp = zeros(P, 3); +Ap = zeros(P, 3); +for p=1:P, + kp = IIRCalcQp(p, P); + [Bp(p,:),Ap(p,:)] = peaking(fs, fc, kp*Q, gain); +endfor + +B = Bp(1,:); +A = Ap(1,:); +if nargout==0 + close + freqz(B,A); + title('Boost Frequency Response') +end + +endfunction + +function [B,A] = peaking(fs, fc, Q, K) + k1=-cos(2*pi*fc/fs); + k2=(1-tan(Q/2))/(1+tan(Q/2)); + Pz = [1 k1*(1+k2) k2]; % define denominator coefficients + Qz = [k2 k1*(1+k2) 1]; % define numerator coefficients + Num = (Pz*(1+K) + Qz*(1-K))/2; + Den = Pz; + B = Num; + A = Den; + +endfunction + +function kp = IIRCalcQp(p, P) + kp = 1.0/(2*sin(pi*(2*p-1)/(2*P))); + +endfunction diff --git a/equalizer/eq_eval.m b/equalizer/eq_eval.m new file mode 100644 index 0000000..8349ca9 --- /dev/null +++ b/equalizer/eq_eval.m @@ -0,0 +1,81 @@ +## Copyright (C) 2020 Jens +## +## This program is free software: you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation, either version 3 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program. If not, see +## . + +## -*- texinfo -*- +## @deftypefn {} {@var{retval} =} eq_eval (@var{input1}, @var{input2}) +## +## @seealso{} +## @end deftypefn + +## Author: Jens +## Created: 2020-06-22 + +function points = eq_eval (filename_in) +P = 1 +fs = 48000 +[x, fs] = audioread(filename_in); +filename_out = to_filename_out(filename_in, "out"); + +points = []; +points = addPoint_peak_notch(points, 120, 60, 6, P); +points = addPoint_peak_notch(points, 540, 480, 3, P); +points = addPoint_peak_notch(points, 1250, 2500, -3, P); +points = addPoint_peak_notch(points, 8000, 8000, 12, P); +points = addPoint_peak_notch(points, 4000, 8000, -12, P); +points = addPoint_peak_notch(points, 15625, 1000, -16, P); + +Np = length(points) +y = x; + +for n = 1:Np, + P = points(n).P; + gain = 10^(points(n).g/20); + + [B,A] = boost3(fs, points(n).f, points(n).b, gain, P); + + for p=1:P, + y = filter(B(p,:), A(p,:), y); + end +endfor +audiowrite(filename_out, y./max(y), fs) + +% Audition +for n = 1:Np, + P = points(n).P; + gain = 10^(points(n).g/20); + + [B,A] = boost3(fs, points(n).f, points(n).b, gain, P); + + for p=1:P, + y = filter(B(p,:), A(p,:), x) - x; + end + suffix = sprintf("%d_%d_%d", points(n).f, points(n).b, points(n).g); + filename_aud = to_filename_out(filename_in, suffix); + + audiowrite(filename_aud, y, fs) + +endfor + +endfunction + +function points_out = addPoint_peak_notch(points_in, fc_Hz, bw_Hz, gain_dB, P) + p.f = fc_Hz; + p.b = bw_Hz; + p.g = gain_dB; + p.P = P; + p.type = "PeakNotch"; + points_out = [points_in p]; +endfunction diff --git a/kalman_eval.m b/kalman_eval.m index 4f2a39a..b5f702f 100644 --- a/kalman_eval.m +++ b/kalman_eval.m @@ -63,7 +63,7 @@ function kalman_eval(varargin) Z = zeros(N,1); % Matrix Z is the measurement noise Y = zeros(N,1); % Matrix Y contains measurement data - H(2,2) = 0; +% H(2,2) = 0; % Known are xp (k|k ), u(k ), P(k|k ) and the new measurement z(k+1). _x = []; diff --git a/stock/diba_aktie.m b/stock/diba_aktie.m new file mode 100644 index 0000000..912251e --- /dev/null +++ b/stock/diba_aktie.m @@ -0,0 +1,10 @@ +function diba_aktie(N, kurs, k_var) + +v_nom = N*kurs +provision_fix = min(69.90, 4.90 + 0.25/100*v_nom) +provision_var = k_var * v_nom + +v_eff = provision_fix + provision_var + v_nom + + +endfunction diff --git a/tempsensor_eval.m b/tempsensor_eval.m new file mode 100644 index 0000000..d81dfdd --- /dev/null +++ b/tempsensor_eval.m @@ -0,0 +1,48 @@ +## Copyright (C) 2019 Jens +## +## This program is free software: you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation, either version 3 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program. If not, see +## . + +## -*- texinfo -*- +## @deftypefn {} {@var{retval} =} tempsensor_eval (@var{input1}, @var{input2}) +## +## @seealso{} +## @end deftypefn + +## Author: Jens +## Created: 2019-04-01 + +function tempsensor_eval () + +data = load('temp.log'); +t = data(:,2); +n = data(:,1); + +dt = [0 diff(t)']'; + +L = 50 +dt2 = []; +for S=1:(length(t)-L) + W = (0:(L-1))'; + t(S:S+(L-1)); + tpf = polyfit(W, t(S:S+(L-1)), 1); + dt2 = [dt2 tpf(1)*60]; +end + +subplot(2,1,1) +plot(n, t); grid; +subplot(2,1,2) +plot(1:length(dt2), dt2); grid; + +endfunction