forked from tszalay/PoreView
-
Notifications
You must be signed in to change notification settings - Fork 0
/
plot_noise.m
84 lines (68 loc) · 2.36 KB
/
plot_noise.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
function plot_noise(sigdata, trange)
%PLOT_NOISE Makes or adds to existing noise plot
% plot_noise(sigdata, trange)
% Takes a SignalData object and a time range as input, uses the same
% algorithm as ClampFit (I think).
% do same thing as ClampFit does - average spectral segs
% ClampFit does 2*65536 pts per seg
% get start and end index
irange = floor(trange/sigdata.si);
fftsize = 2*65536;
% only process real signals
dfft = zeros(fftsize,sigdata.nsigs);
% number of frames
nframes = 0;
wh = waitbar(0,'Calculating power spectrum...','Name','CrampFit');
for ind=irange(1):fftsize:irange(2)
% get only the real signals
d = sigdata.get(ind:ind+fftsize-1,1+(1:sigdata.nsigs));
% quit if we don't have enough points
if size(d,1) < fftsize
break
end
% calculate power spectrum
df = sigdata.si*abs(fft(d)).^2/fftsize;
% and add to fft accum.
dfft = dfft + df;
nframes = nframes + 1;
waitbar((ind-irange(1))/(irange(2)-irange(1)));
end
close(wh);
% do the averaging
dfft = dfft / nframes;
% and calculate the frequency range
freqs = 1/sigdata.si*(0:fftsize-1)/fftsize;
% range of freqs to plot, only do half (Nyquist and all)
imax = floor(size(dfft,1)/2);
hf = findobj('Name','Noise Power Spectrum');
if isempty(hf)
% make a new plot figure, leave menu bar
hf = figure('Name','Noise Power Spectrum','NumberTitle','off');
% Build a new color map
CO = [ 0.0 0.2 0.7;...
0.0 0.5 0.0;...
0.7 0.1 0.1;...
0.0 0.6 0.6;...
0.5 0.0 0.5;...
0.6 0.6 0.3 ];
% Set the color order
set(hf, 'DefaultAxesColorOrder', CO);
% and make axes
ax = axes('Parent',hf,'XScale','log','YScale','log',...
'XLimMode','manual','YLimMode','auto','NextPlot','add',...
'TickDir','out');
% scale x axis
set(ax,'XLim',[1 freqs(imax)]);
% grid
grid on
grid minor
% labels
title('Noise Power Spectrum');
ylabel('Power (nA^2/Hz)')
xlabel('Frequency (Hz)')
end
% bring figure to front
figure(hf);
% and plot
plot(freqs(1:imax)',dfft(1:imax,:));
end