jens/matlab
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

- added

Browse Source 
git-svn-id: http://moon:8086/svn/matlab/trunk@128 801c6759-fa7c-4059-a304-17956f83a07c
This commit is contained in:
jens
2020-07-23 11:10:12 +00:00
parent 0631c89f10
commit 1f5b0b0570
8 changed files with 288 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files
equalizer/boost.m
+28
View File
@@ -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
equalizer/boost2.m
+29
View File
@@ -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 <rbj@audioimagination.com>
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
equalizer/boost3.m
+48
View File
@@ -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 <rbj@audioimagination.com>
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
equalizer/boost3a.m
+43
View File
@@ -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 <rbj@audioimagination.com>
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
equalizer/eq_eval.m
+81
View File
@@ -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
## <https://www.gnu.org/licenses/>.
## -*- texinfo -*-
## @deftypefn {} {@var{retval} =} eq_eval (@var{input1}, @var{input2})
##
## @seealso{}
## @end deftypefn
## Author: Jens <jens@orion>
## 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
kalman_eval.m
+1 -1
View File
@@ -63,7 +63,7 @@ function kalman_eval(varargin)
Z = zeros(N,1); % Matrix Z is the measurement noise Z = zeros(N,1); % Matrix Z is the measurement noise
Y = zeros(N,1); % Matrix Y contains measurement data 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). % Known are xp (k|k ), u(k ), P(k|k ) and the new measurement z(k+1).
_x = []; _x = [];
stock/diba_aktie.m
+10
View File
@@ -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
tempsensor_eval.m
+48
View File
@@ -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
## <https://www.gnu.org/licenses/>.
## -*- texinfo -*-
## @deftypefn {} {@var{retval} =} tempsensor_eval (@var{input1}, @var{input2})
##
## @seealso{}
## @end deftypefn
## Author: Jens <jens@orion>
## 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