-
Notifications
You must be signed in to change notification settings - Fork 0
/
findbodedata2.m
81 lines (72 loc) · 3.69 KB
/
findbodedata2.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
function [gain, phase, w] = findbodedata2(input, output, Ts, neglect)
% function [gain, phase, w] = findbodedata2(input, output, Ts, neglect)
% This function receives a response (output) to a sinusoid (input),
% both sampled with sampling time Ts (in s, default is 1 s),
% and calculates gain (dB), phase lag (º), and input and output frequency (w, in rad/s).
% The first neglect(1) points (default is 0) and the last neglect(2) points (default is 0) are not considered.
% Duarte Valério 2007-2013; see below for code incorporated from Matlab FileExchange
if nargin < 4 % by default no points are neglected anywhere
neglect = [0 0];
end
if length(neglect) == 1 % by default no points are neglected at the end
neglect = [neglect 0];
end
if length(input) < sum(neglect)
error('You cannot neglect all the data.')
end
if length(input) ~= length(output)
error('The two sinusoids must have the same lengths.')
end
if nargin < 3
Ts = 1;
end
input = input(neglect(1)+1:end-neglect(2));
output = output(neglect(1)+1:end-neglect(2));
input = input - mean(input);
output = output - mean(output);
gain = 20*log10( sqrt(mean(output.^2)) / sqrt(mean(input.^2)) );
lengthSignal = length(input); % number of samples
windowsize = 2^nextpow2(lengthSignal);
INPUT = fft(input, windowsize) / lengthSignal; % Fourier Transform of input
[~, indexINPUT] = max( abs(INPUT(1 : windowsize/2+1)) );
phaseINPUT = angle(INPUT(indexINPUT)); % phase of the most significant frequency component of input
OUTPUT = fft(output, windowsize) / lengthSignal; % Fourier Transform of output
[~, indexOUTPUT] = max( abs(OUTPUT(1 : windowsize/2+1)) );
phaseOUTPUT = angle(OUTPUT(indexOUTPUT)); % phase of the most significant frequency component of output
phase = rad2deg(phaseOUTPUT - phaseINPUT);
f = (2*pi/Ts) * (0:windowsize/2) / windowsize;
w = [f(indexINPUT), f(indexOUTPUT)];
% ORIGINAL HEADER AND COPYRIGHT OF FILE GIVEPD.M FROM MATLAB FILE EXCHANGE, ADAPTED INTO LINES 34-41
% Filename: givepd.m
% Author: Shashank G. Sawant
% email: sgsawant@gmail.com
% Description: Given 2 sinusoidal signals of the
% same frequency, the function gives the "phase difference" between the
% 2 given signals
% The phase difference is in RADIANS!!!!!!
% The output is limited to pd={-pi,pi}radians
% Time Stamp: 2010 October 19 2043hrs
% Copyright (c) 2010, Shashank Sawant
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
% * Redistributions of source code must retain the above copyright
% notice, this list of conditions and the following disclaimer.
% * Redistributions in binary form must reproduce the above copyright
% notice, this list of conditions and the following disclaimer in
% the documentation and/or other materials provided with the distribution
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.