-
Notifications
You must be signed in to change notification settings - Fork 2
/
polariz.m
executable file
·125 lines (113 loc) · 4.22 KB
/
polariz.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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% PROGRAM:
% polariz.m
%
% PROGRAMMER:
% Matt Haney
%
% Last revision date:
% 22 May 2009
% Last modified by Dylan Mikesell (19 December 2016)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Program polariz is a Matlab script that compares polarization
% estimates based on covariance and coherence matrices
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this is a script
clear all
close all
clc
addpath('./functions');
%% Step 1: Make synthetic 3-component times series
%--------------------------------------------------------------------------
% length of synthetic time series
tln = 10000;
%--------------------------------------------------------------------------
% dominant frequency of synthetic time series in Hz
f = 2; % [Hz]
omga = 2*pi*f; % [radians/s]
%--------------------------------------------------------------------------
% time sample interval
dt = 0.001; % [s]
%--------------------------------------------------------------------------
% Determine window to calculate polarization over, in samples
cycs = 2; % number of cycles (2 to 3 is usually sufficient)
wndo = floor( (1/f) * (1/dt) ) * cycs; % samples per cycle times # of cycles
%--------------------------------------------------------------------------
[dtac,tt] = makeSynthetic(tln,dt,omga); % compute the synthetics
%--------------------------------------------------------------------------
% plot the three-component synthetic data
lSize = 2;
h = figure('Color','white');
subplot(3,1,1)
plot(tt,dtac(1,:),'LineWidth',lSize); axis([tt(1) tt(end) -3.5 3.5]); grid on;
title('Three-component synthetic data');
ylabel('Vertical [Z]');
subplot(3,1,2)
plot(tt,dtac(2,:),'LineWidth',lSize); axis([tt(1) tt(end) -3.5 3.5]); grid on;
ylabel('East [E]');
subplot(3,1,3)
plot(tt,dtac(3,:),'LineWidth',lSize); axis([tt(1) tt(end) -3.5 3.5]); grid on;
ylabel('North [N]'); xlabel('Time [s]');
%--------------------------------------------------------------------------
% Fix figure properties
fsize = 16;
% set( findall( h, '-property', 'FontWeight' ), 'FontWeight', 'Bold' );
set( findall( h, '-property', 'Fontsize' ), 'Fontsize', fsize );
%
%% Apply polar coherency method
%
% calculate the polarization parameters using the complex-valued coherency
% method
[azim, incd, ellip] = polar_coherency( dtac, wndo );
%--------------------------------------------------------------------------
% plot the coherency-based polarization parameters
lSize = 2;
h = figure('Color','white');
subplot(3,1,1);
plot(tt,azim,'LineWidth',lSize); %axis([tt(1) tt(end) 60 75]);
grid on;
axis([tt(1) tt(end) 50 100]);
title('Coherency matrix method (Vidale, BSSA, 1986)');
ylabel('Azimuth [deg]');
subplot(3,1,2)
plot(tt,incd,'LineWidth',lSize); axis([tt(1) tt(end) 50 100]); grid on;
ylabel('Incididence [deg]');
subplot(3,1,3)
plot(tt,ellip,'LineWidth',lSize); axis([tt(1) tt(end) 0 1]); grid on;
xlabel('Time [s]'); ylabel('Ellipticity');
%--------------------------------------------------------------------------
% Fix figure properties
fsize = 16;
% set( findall( h, '-property', 'FontWeight' ), 'FontWeight', 'Bold' );
set( findall( h, '-property', 'Fontsize' ), 'Fontsize', fsize );
%
%% Apply polar covariance method
%
% calculate the polarization parameters using the real-valued covariance
% method
[azim, incd, ellip] = polar_covariance( dtac, wndo );
%--------------------------------------------------------------------------
% plot the covariance-based polarization parameters
lSize = 2;
h = figure('Color','white');
subplot(3,1,1)
plot(tt,azim,'LineWidth',lSize); %axis([tt(1) tt(end) 60 75]);
grid on;
axis([tt(1) tt(end) 50 100]);
title('Covariance matrix method');
ylabel('Azimuth [deg]');
subplot(3,1,2)
plot(tt,incd,'LineWidth',lSize); axis([tt(1) tt(end) 50 100]); grid on;
ylabel('Incididence [deg]');
subplot(3,1,3)
plot(tt,ellip,'LineWidth',lSize); axis([tt(1) tt(end) 0 1]); grid on;
xlabel('Time [s]'); ylabel('Ellipticity');
%--------------------------------------------------------------------------
% Fix figure properties
fsize = 16;
% set( findall( h, '-property', 'FontWeight' ), 'FontWeight', 'Bold' );
set( findall( h, '-property', 'Fontsize' ), 'Fontsize', fsize );