1function [R, Rax, F] = gradSpectrum(obj, FB, fmax, doPlot)
2% function [R, Rax, F] = seq.gradSpectrum(FB, [fmax = 3000], [doPlot = true])
3%
4% Get (and optionally plot) frequency response of a Pulseq sequence.
5%
6% Input 'FB' is either a struct array with the forbidden bands,
7% or the name of a Siemens ASC file.
8%
9% Inputs
10% seq Pulseq sequence object
11% FB [num_bands] Forbidden frequency bands (struct array)
12% FB(1).freq center frequency of first forbidden band
13% FB(1).bw bandwidth of first forbidden band
14% FB string Siemens ASC file name
15% fmax [1] Max frequency range
16% doPlot TRUE/false Plot or just return values
17%
18% Outputs
19% R [n] frequency response (root-sum-of-squares of all axes)
20% Rax [3 n] frequency response for individual axes
21% F [n] frequency locations (Hz)
22%
23% Function version of demoUnsorted/gradSpectrum.m
25% defaults
26if nargin < 4
27 doPlot = true;
28end
29if nargin < 3
30 fmax = 3000;
31end
32if ~exist('FB','var')
33 FB=[];
34end
36% Read ASC file if provided
37if ischar(FB)
38 ascData=mr.Siemens.readasc(FB);
39 clear FB
40 for i=1:length(ascData.asGPAParameters(1).sGCParameters.aflAcousticResonanceFrequency)
41 if (ascData.asGPAParameters(1).sGCParameters.aflAcousticResonanceFrequency(i)>0)
42 FB(i).freq = ascData.asGPAParameters(1).sGCParameters.aflAcousticResonanceFrequency(i);
43 FB(i).bw = ascData.asGPAParameters(1).sGCParameters.aflAcousticResonanceBandwidth(i);
44 end
45 end
46end
48% Calculate spectrum/spectrogramm
49dt=obj.sys.gradRasterTime; % time raster
50nwin=5000; % 0.05s
51os=3; % frequency oversampling for prettier peaks
53faxis=(0:(nwin/2-1))/nwin/dt/os;
54nfmax=sum(faxis<=fmax);
56wave_data=obj.waveforms_and_times();
57ng=length(wave_data);
58tmax=0;
59for i=1:ng
60 if ~isempty(wave_data{i})
61 tmax=max(tmax, wave_data{i}(1,end));
62 end
63end
64if tmax==0
65 error('Empty sequence passed to gradSpectrum()');
66end
67nt=ceil(tmax/dt);
68tmax=nt*dt;
70gw=zeros(ng,nt);
71for i=1:ng
72 gw(i,:)=interp1(wave_data{i}(1,:),wave_data{i}(2,:),((1:nt)-0.5)*dt,'linear',0);
73 % alternative (to be checked in the future)
74 % it is actually much more appropriate to calculate the spectrium of
75 % the derivative(!) of the gradient wave form and not the waveform
76 % itself, at least for the cound/noise of the gradients...
77 %gw(i,1:end-1)=diff(interp1(wave_data{i}(1,:),wave_data{i}(2,:),((1:nt)-0.5)*dt,'linear',0));
78end
80gs=[];
82ng=size(gw,1);
83for g=1:ng
84 x=gw(g,:);
85 nx = length(x);
87 nx=ceil(nx/nwin)*nwin;
88 if nx>length(x)
89 x=[x, zeros(1,nx-length(x))]; % zerofill
90 end
92 nseg1=nx/nwin;
93 xseg=zeros(nseg1*2-1,nwin*os);
95 xseg(1:2:end,1:nwin)=reshape(x,[nwin,nseg1])';
96 if nseg1>1
97 xseg(2:2:end,1:nwin)=reshape(x(1+nwin/2:end-nwin/2),[nwin,nseg1-1])';
98 end
100 xseg_dc=mean(xseg,2);
101 xseg=xseg-xseg_dc(:,ones(1,nwin*os));
103 if nseg1>1 % WARNING: this introduces an inconsistency between short and long sequences in term os the peak amplitudes
104 cwin=0.5*(1-cos(2*pi*(1:nwin)/nwin));
105 xseg(:,1:nwin)=xseg(:,1:nwin).*cwin(ones(size(xseg,1),1),:);
106 end
108 fseg=abs(fft(xseg,[],2));
109 fseg=fseg(:,1:end/2);
111 if nseg1>1
112 gs = [gs; mean(fseg.^2).^0.5]; % sos
113 %figure; plot(faxis(1:nfmax),sum(fseg(:,1:nfmax).^2).^0.5);
114 else
115 gs = [gs; abs(fseg)]; % add abs
116 end
117end
119% Define return values
120F = faxis(1:nfmax);
121Rax = gs(:,1:nfmax);
122R = sum(gs(:,1:nfmax).^2).^0.5;
124% plot
125if ~doPlot
126 return;
127end
129figure; plot(faxis(1:nfmax),gs(:,1:nfmax));
130hold on; plot(faxis(1:nfmax),sum(gs(:,1:nfmax).^2).^0.5); % sos
131xlabel('frequency / Hz');
133% alternative "stained glass" plots
134% clr = repmat('rgb', [1 5]);
135% if ~isempty(FB)
136% for i=1:length(FB)
137% if FB(i).freq > 0
138% l = FB(i).freq-FB(i).bw/2;
139% r = FB(i).freq+FB(i).bw/2;
140% t = max(R);
141% b = 0;
142% h = fill([l r r l], [b b t t], clr(i));
143% h.FaceAlpha = 0.15;
144% end
145% end
146% end
148if ~isempty(FB)
149 for i=1:length(FB)
150 if FB(i).freq > 0
151 if exist('xline','file')
152 xline(FB(i).freq,'-');
153 xline(FB(i).freq-FB(i).bw/2,'--');
154 xline(FB(i).freq+FB(i).bw/2,'--');
155 else
156 % Octave fallback
157 yl=ylim;
158 line([1 1]*FB(i).freq, yl, 'Color','k','LineStyle','-');
159 line([1 1]*(FB(i).freq-FB(i).bw/2), yl, 'Color','k','LineStyle','--');
160 line([1 1]*(FB(i).freq+FB(i).bw/2), yl, 'Color','k','LineStyle','--');
161 end
162 end
163 end
164end
166legend({'Gx','Gy','Gz','Gtot'});
168%nov = floor(nsc/2);
169%nff = max(256,2^nextpow2(nsc));
171%t = spectrogram(x,hamming(nsc),nov);%,nff);
172%t = spectrogram(x,rectwin(nsc),nov);
173%maxerr = max(abs(abs(t(:))-abs(s(:))))