1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
|
function [ b, a ] = pei_tseng_notch ( frequencies, bandwidths )
//Return coefficients for an IIR notch-filter.
//Calling Sequence:
//[b, a] = pei_tseng_notch(frequencies, bandwidths)
//Parameters:
//frequencies: Real scalar/vector representing filter frequencies.
//bandwidths: Real scalar scalar/vector representing bandwidths to be used with filter.
//Description:
//THis function returns coefficients for an IIR notch-filter with one or more filter frequencies and according bandwidths.
//The filter is based on a all pass filter that performs phasereversal at filter frequencies. This leads to removal of those frequencies of the original and phase-distorted signal.
//Examples:
//sf = 800; sf2 = sf/2;
//data = [[1;zeros(sf-1,1)],sinetone(49,sf,1,1),sinetone(50,sf,1,1),sinetone(51,sf,1,1)];
//[b,a] = pei_tseng_notch ( 50 / sf2, 2/sf2 )
//b = 0.99213 -1.83322 0.99213
//a = 1.00000 -1.83322 0.98426
if (nargin() ~= 2)
error("Wrong number of input arguments.");
elseif ( ~isvector(frequencies) | ~isvector(bandwidths) )
if ( ~isscalar(frequencies) | ~isscalar(bandwidths) )
error("All arguments must be vectors or scalars.");
end
elseif ( length(frequencies) ~= length (bandwidths) )
error("All arguments must be of equal length.");
elseif (~and( frequencies > 0 & frequencies < 1 ))
error("All frequencies must be in (0, 1).");
elseif (~and( bandwidths > 0 & bandwidths < 1 ))
error("All bandwidths must be in (0, 1).");
end
frequencies = frequencies (:)';
bandwidths = bandwidths (:)';
frequencies = frequencies*%pi;
bandwidths = bandwidths*%pi;
M2 = 2 * length(frequencies);
a = [ frequencies - bandwidths / 2; frequencies ];
omega = a(:);
factors = ( 1 : 2 : M2 );
b = [ -%pi * factors + %pi / 2; -%pi * factors ];
phi = b(:);
t_beta = tan ( ( phi + M2 * omega ) / 2 );
Q = zeros (M2, M2);
for k = 1 : M2
Q ( :, k ) = sin ( k .* omega ) - t_beta .* cos ( k .* omega );
end
h_a = ( Q \ t_beta )';
denom = [ 1, h_a ];
num = [ flipdim(h_a, 2), 1 ];
a = denom;
b = ( num + denom ) / 2;
endfunction
//input validation:
//assert_checkerror("pei_tseng_notch()", "Wrong number of input arguments.");
//assert_checkerror("pei_tseng_notch([1, 2, 3])", "Wrong number of input arguments.");
//assert_checkerror("pei_tseng_notch([1, 2; 3, 4], [4; 5; 6])", "All arguments must be vectors or scalars.");
//assert_checkerror("pei_tseng_notch([1, 2, 3, 4], [5, 6])", "All arguments must be of equal length.");
//assert_checkerror("pei_tseng_notch([1, 2, 3], [4, 5, 6])", "All frequencies must be in (0, 1).");
//assert_checkerror("pei_tseng_notch([0.1, 0.2, 0.3], [4, 5, 6])", "All bandwidths must be in (0, 1).");
//tests:
//assert_checkequal(pei_tseng_notch(0.2, 0.4), [0 0 0]);
//assert_checkalmostequal(pei_tseng_notch([0.2, 0.5, 0.9], [0.1, 0.3, 0.5]), [0.41549, 0.11803, -0.03227, 0.23606, -0.03227, 0.11803, 0.41549], 5*10^-4);
//assert_checkequal(pei_tseng_notch([0.2; 0.5; 0.9], [0.1, 0.3, 0.5]), pei_tseng_notch([0.2, 0.5, 0.9], [0.1; 0.3; 0.5]));
//assert_checkalmostequal(pei_tseng_notch([0.2, 0.7], [0.1, 0.3]), [0.60331, -0.26694, 0.05905, -0.26694, 0.60331], 5*10^-4);
|
