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
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
|
/*Calling Sequence
h = hilbert1 (f)
h = hilbert1 (f, N)
h = hilbert1 (f, N, dim)
Description
Analytic extension of real valued signal.
h = hilbert (f) computes the extension of the real valued signal f to an analytic signal.
If f is a matrix, the transformation is applied to each column.
For N-D arrays, the transformation is applied to the first non-singleton dimension.
real (h) contains the original signal f. imag (h) contains the Hilbert transform of f.
hilbert (f, N) does the same using a length N Hilbert transform. The result will also have length N.
hilbert (f, [], dim) or hilbert (f, N, dim) does the same along dimension dim.
Dependencies
fft1, ifft1, ipermute
Example
//the magnitude of the hilbert transform eliminates the carrier
t=linspace(0,10,1024);
x=5*cos(0.2*t).*sin(100*t);
plot(t,x,t,abs(hilbert(x)));
*/
function f=hilbert1(f, N ,dim )
// ------ PRE: initialization and dimension shifting ---------
nargin = argn(2);
if (nargin<1 || nargin>3)
error("Please enter valid number of inputs")
end
if ~isreal(f)
warning ('HILBERT: ignoring imaginary part of signal');
f = real (f);
end
D=ndims(f);
select nargin
case 1 then
N=[];
dim=[];
case 2 then
dim=[]
end
// Dummy assignment.
order=1;
if isempty(dim)
dim=1;
if sum(size(f)>1)==1
// We have a vector, find the dimension where it lives.
dim=find(size(f)>1);
end
else
if (length(dim)~=1 || ~or(type(dim)==[1 5 8]))
error('HILBERT: dim must be a scalar.');
end
if modulo(dim,1)~=0
error('HILBERT: dim must be an integer.');
end
if (dim<1) || (dim>D)
error('HILBERT: dim must be in the range from 1 to %d.',D);
end
end
if (length(N)>1 || ~or(type(N)==[1 5 8]))
error('N must be a scalar.');
elseif (~isempty(N) && modulo(N,1)~=0)
error('N must be an integer.');
end
if dim>1
order=[dim, 1:dim-1,dim+1:D];
// Put the desired dimension first.
f=permute(f,order);
end
Ls=size(f,1);
// If N is empty it is set to be the length of the transform.
if isempty(N)
N=Ls;
end
// moduloember the exact size for later and modify it for the new length
permutedsize=size(f);
permutedsize(1)=N;
// Reshape f to a matrix.
f=resize_matrix(f,size(f,1),length(f)/size(f,1));
W=size(f,2);
if ~isempty(N)
siz=size(f);
siz(1)=N;
f=resize_matrix(f,siz);
end
// ------- actual computation -----------------
if N>2
f=fft1(f);
if modulo(N,2)==0
f=[f(1,:);
2*f(2:N/2,:);
f(N/2+1,:);
zeros(N/2-1,W)];
else
f=[f(1,:);
2*f(2:(N+1)/2,:);
zeros((N-1)/2,W)];
end
f=ifft1(f);
end
// ------- POST: Restoration of dimensions ------------
// Restore the original, permuted shape.
f=matrix(f,permutedsize);
if dim>1
// Undo the permutation.
f=ipermute(f,order);
end
endfunction
|
