-
Notifications
You must be signed in to change notification settings - Fork 6
/
demo1.m
145 lines (97 loc) · 2.69 KB
/
demo1.m
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
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
clear all
close all
clc
% create time series grid as graph
n=101;
W=zeros(n,n);
for i=1:n-1
W(i,i+1)=1;
W(i+1,i)=1;
end
W(1,n)=1;
W(n,1)=1;
% calculate combinatorial Laplacian Matrix
d = sum(W,2);
L = diag(d)-W;
% calculate basis
[u v]=eig(L);
% s=sign(u(1,:));
% s(s==0)=1;
% u=u.*s;
figure;subplot(1,3,1);plot(u(:,1));ylim([-0.15 0.15]);xlim([1 101]);
title(['1st eigenvector, value=' num2str(v(1,1))]);
subplot(1,3,2);plot(u(:,2));ylim([-0.15 0.15]);xlim([1 101]);
title(['2nd eigenvector, value=' num2str(v(2,2))]);
subplot(1,3,3);plot(u(:,3));ylim([-0.15 0.15]);xlim([1 101]);
title(['3th eigenvector, value=' num2str(v(3,3))]);
% make eignevalue as vector
v=diag(v);
% get maximum eigenvalue
lmax=max(v);
% create arbitrary signal
s=randn(size(W,1),1);
%s(50)=1;
% calcualte graph frequency profile
f=u'*s;
% calculate classical frequency profile
cf=abs(fft(s));
cf=cf(1:51);
% F=[0;f];
% F=reshape(F,2,51);
% f=sum(abs(F));
figure;subplot(1,2,1);plot(v,f,'r*-');%plot(v(1:2:end),f,'r*-');
xlabel('eigenvalues');ylabel('coefficient');
title('Graph Frequency profile');
subplot(1,2,2);plot(0:50,cf,'b*-');
xlabel('Frequencies');ylabel('magnitude');
title('Frequency profile');
% create filter for classical signal processing
fs=fftshift(fft(s));
f=linspace(-(n-1)/2,(n-1)/2,n)';
flt = exp(-abs(f)*0.1);
% apply that filter on to time series signal
ot=ifft(ifftshift(flt.*fs));
% create the same filter for graph signal processing
flt=[1:(n-1)/2]';
flt=[flt flt]';
flt=[0; flt(:)];
flt = exp(-abs(flt)*0.1);
% apply that signal on to graph signal
sf=u*(flt.*(u'*s));
figure;
hold on;plot(s,'k--');plot(ot,'b-','linewidth',2);xlim([1 n]);
plot(sf,'r-','linewidth',1);xlabel('node order or time');
plot(ot-sf,'g-')
legend({'original signal','graph filtering','time series filtering','differences'});
%% show graph and signal
run gspbox/gsp_start
coor=u(:,2:3);
G=gsp_graph(W,coor);
figure;gsp_plot_signal(G,s)
title('Input signal');
figure;gsp_plot_signal(G,sf)
title('Filtered signal');
figure;
frm=0;
for i=1:size(u,2)
frm=frm+1;
ff=abs(fft(u(:,i)));
[a b]=max(ff(1:(length(ff)+1)/2));
msg=[num2str(i) 'th eigenvector measured freq:' num2str(b-1) 'Hz'];
h=figure(5);
plot(u(:,i));
xlabel('sample suppose total 1 sec');
ylim([-0.2 0.2]);
xlim([1 n])
title(msg);
%pause(0.1);
frame = getframe(h);
im = frame2im(frame);
[imind,cm] = rgb2ind(im,256);
%Write to the GIF File
if frm == 1
imwrite(imind,cm,'eigenvectors11.gif','gif', 'Loopcount',inf);
else
imwrite(imind,cm,'eigenvectors11.gif','gif','WriteMode','append');
end
end