Category: Matlab
Category Archives: Matlab
I don’t know what I’m wrong..
years = [1000, 1650, 1800, 1900, 1950, 1960, 1970, 1980, 1990];
population = [0.34, 0.545, 0.907, 1.61, 2.51, 3.15, 3.65, 4.2, 5.3];
x = years;
y = population;
A = [ones(length(x), 1), x’];
b = y’ ./ (1 + x’);
[Q, R] = householderQR(A);
c1 = Q’ * b;
c1 = c1(1:2);
x_hh = R c1;
tt = linspace(min(x), max(x), 100);
yy = (x_hh(1) + x_hh(2)*tt) ./ (1 + tt);
clf;
hold on;
plot(tt, yy, ‘LineWidth’, 2);
plot(x, y, ‘k*’, ‘MarkerSize’, 10, ‘LineWidth’, 2);
hold off;
grid on;
legend(‘Fitting Curve’, ‘Data Points’);
title(‘World Population Fitting Using Householder QR Factorization’);
xlabel(‘Year’);
ylabel(‘Population (billions)’);
function [Q, R] = householderQR(A)
[m, n] = size(A);
Q = eye(m);
R = A;
for k = 1:n
% Step 1: Define alpha
alpha = -sign(R(k, k)) * norm(R(k:m, k));
% Step 2: Define v
v = R(k:m, k);
v(1) = v(1) – alpha;
v = v / norm(v);
% Step 3: Define H
H = eye(m-k+1) – 2 * (v * v’);
H_full = eye(m);
H_full(k:m, k:m) = H;
% Step 4: Update R and Q
R = H_full * R;
Q = Q * H_full’;
end
% Adjust R to be upper triangular of size n x n
R = R(1:n, :);
end
What I wrote is the code of the below question.
A quick way of finding a function of the form f(x) ≈ (a + bx) / (1 + cx) is to apply the least squares method to the problem (1+cx)f(x) ≈ (a+bx). Use this technique to fit the world population (billions) data given using Householder QR method (without using built-in functions). Make a plot of the original data points along with resulting curve.
I don’t know why my code of my plot is not working. Could someone modify my code to work well?years = [1000, 1650, 1800, 1900, 1950, 1960, 1970, 1980, 1990];
population = [0.34, 0.545, 0.907, 1.61, 2.51, 3.15, 3.65, 4.2, 5.3];
x = years;
y = population;
A = [ones(length(x), 1), x’];
b = y’ ./ (1 + x’);
[Q, R] = householderQR(A);
c1 = Q’ * b;
c1 = c1(1:2);
x_hh = R c1;
tt = linspace(min(x), max(x), 100);
yy = (x_hh(1) + x_hh(2)*tt) ./ (1 + tt);
clf;
hold on;
plot(tt, yy, ‘LineWidth’, 2);
plot(x, y, ‘k*’, ‘MarkerSize’, 10, ‘LineWidth’, 2);
hold off;
grid on;
legend(‘Fitting Curve’, ‘Data Points’);
title(‘World Population Fitting Using Householder QR Factorization’);
xlabel(‘Year’);
ylabel(‘Population (billions)’);
function [Q, R] = householderQR(A)
[m, n] = size(A);
Q = eye(m);
R = A;
for k = 1:n
% Step 1: Define alpha
alpha = -sign(R(k, k)) * norm(R(k:m, k));
% Step 2: Define v
v = R(k:m, k);
v(1) = v(1) – alpha;
v = v / norm(v);
% Step 3: Define H
H = eye(m-k+1) – 2 * (v * v’);
H_full = eye(m);
H_full(k:m, k:m) = H;
% Step 4: Update R and Q
R = H_full * R;
Q = Q * H_full’;
end
% Adjust R to be upper triangular of size n x n
R = R(1:n, :);
end
What I wrote is the code of the below question.
A quick way of finding a function of the form f(x) ≈ (a + bx) / (1 + cx) is to apply the least squares method to the problem (1+cx)f(x) ≈ (a+bx). Use this technique to fit the world population (billions) data given using Householder QR method (without using built-in functions). Make a plot of the original data points along with resulting curve.
I don’t know why my code of my plot is not working. Could someone modify my code to work well? years = [1000, 1650, 1800, 1900, 1950, 1960, 1970, 1980, 1990];
population = [0.34, 0.545, 0.907, 1.61, 2.51, 3.15, 3.65, 4.2, 5.3];
x = years;
y = population;
A = [ones(length(x), 1), x’];
b = y’ ./ (1 + x’);
[Q, R] = householderQR(A);
c1 = Q’ * b;
c1 = c1(1:2);
x_hh = R c1;
tt = linspace(min(x), max(x), 100);
yy = (x_hh(1) + x_hh(2)*tt) ./ (1 + tt);
clf;
hold on;
plot(tt, yy, ‘LineWidth’, 2);
plot(x, y, ‘k*’, ‘MarkerSize’, 10, ‘LineWidth’, 2);
hold off;
grid on;
legend(‘Fitting Curve’, ‘Data Points’);
title(‘World Population Fitting Using Householder QR Factorization’);
xlabel(‘Year’);
ylabel(‘Population (billions)’);
function [Q, R] = householderQR(A)
[m, n] = size(A);
Q = eye(m);
R = A;
for k = 1:n
% Step 1: Define alpha
alpha = -sign(R(k, k)) * norm(R(k:m, k));
% Step 2: Define v
v = R(k:m, k);
v(1) = v(1) – alpha;
v = v / norm(v);
% Step 3: Define H
H = eye(m-k+1) – 2 * (v * v’);
H_full = eye(m);
H_full(k:m, k:m) = H;
% Step 4: Update R and Q
R = H_full * R;
Q = Q * H_full’;
end
% Adjust R to be upper triangular of size n x n
R = R(1:n, :);
end
What I wrote is the code of the below question.
A quick way of finding a function of the form f(x) ≈ (a + bx) / (1 + cx) is to apply the least squares method to the problem (1+cx)f(x) ≈ (a+bx). Use this technique to fit the world population (billions) data given using Householder QR method (without using built-in functions). Make a plot of the original data points along with resulting curve.
I don’t know why my code of my plot is not working. Could someone modify my code to work well? householderqr, least square method MATLAB Answers — New Questions
help changing subplot sequence in rdmseed
How can i change the subplots sequence in rdmseed.m, while plotting a miniseed file with e.g. 5 stations?
function varargout = rdmseed(varargin)
if nargin > 6
error(‘Too many input arguments.’)
end
% global variables shared with sub-functions
global f fid offset le ef wo rl forcebe verbose notc force
% default input arguments
makeplot = 0; % make plot flag
verbose = 0; % verbose flag/level
forcebe = 0; % force big-endian
ef = 10; % encoding format default
wo = 1; % word order default
rl = 2^12; % record length default
force = 0; % force input argument over blockette 1000 (UNDOCUMENTED)
notc = 0; % force no time correction (over ActivityFlags)
nullhead = 0; % allow null bytes before header
if nargin < 1
[filename,pathname] = uigetfile(‘*’,’Please select a miniSEED file…’);
f = fullfile(pathname,filename);
else
f = varargin{1};
end
if ~ischar(f) || ~exist(f,’file’)
error(‘File %s does not exist.’,f);
end
if nargin > 1
verbose = any(strcmpi(varargin,’v’)) + 2*any(strcmpi(varargin,’vv’)) …
+ 3*any(strcmpi(varargin,’vvv’));
makeplot = any(strcmpi(varargin,’plot’));
forcebe = any(strcmpi(varargin,’be’));
notc = any(strcmpi(varargin,’notc’));
force = any(strcmpi(varargin,’force’));
nullhead = any(strcmpi(varargin,’nullhead’));
end
nargs = (makeplot>0) + (verbose>0) + (forcebe>0) + (notc>0) + (force>0) …
+ (nullhead>0);
if nargin > (1 + nargs)
ef = varargin{2};
if ~isnumeric(ef) || ~any(ef==[0:5,10:19,30:33])
error(‘Argument ENCODINGFORMAT must be a valid FDSN code value.’);
end
end
if nargin > (2 + nargs)
wo = varargin{3};
if ~isnumeric(wo) || (wo ~= 0 && wo ~= 1)
error(‘Argument WORDORDER must be 0 or 1.’);
end
end
if nargin > (3 + nargs)
rl = varargin{4};
if ~isnumeric(rl) || rl < 256 || rem(log(rl)/log(2),1) ~= 0
error(‘Argument RECORDLENGTH must be a power of 2 and greater or equal to 256.’);
end
end
if nargout == 0
makeplot = 1;
end
% sensible limits for multiplexed files
max_channels = 20; % absolute max number of channels to plot
max_channel_label = 6; % max. number of channels for y-labels
% file is opened in Big-Endian encoding (this is encouraged by SEED)
fid = fopen(f,’rb’,’ieee-be’);
le = 0;
offset = 0;
% — tests if the header is mini-SEED
% the 7th character must be one of the "data header/quality indicator", usually ‘D’
header = fread(fid,20,’*char’);
if ~ismember(header(7),’DRMQ’)
if ismember(header(7),’VAST’)
error(‘File seems to be a SEED Volume. Cannot read it.’);
else
if header(1)==0
if nullhead
if verbose
fprintf(‘Null header option: bypassing…’);
end
c = 0;
fseek(fid,0,’bof’);
while c==0
c = fread(fid,1,’*char’);
offset = offset + 1;
end
if verbose
fprintf(‘ %d null bytes.n’,offset);
end
header = fread(fid,6,’*char’);
if ~ismember(header(6),’DRMQ’)
error(‘File is not in mini-SEED format. Cannot read it.’);
else
offset = offset – 1;
end
else
error(‘File starts with null bytes… if you believe it is still a miniseed file, try the ”nullhead” option.’);
end
else
error(‘File is not in mini-SEED format. Cannot read it.’);
end
end
end
i = 1;
while offset >= 0
X(i) = read_data_record;
i = i + 1;
end
fclose(fid);
if nargout > 0
varargout{1} = X;
end
% — analyses data
if makeplot || nargout > 1
% test if the file is multiplexed or a single channel
un = unique(cellstr(char(X.ChannelFullName)));
nc = numel(un);
for i = 1:nc
k = find(strcmp(cellstr(char(X.ChannelFullName)),un{i}));
I(i).ChannelFullName = X(k(1)).ChannelFullName;
I(i).XBlockIndex = k;
I(i).ClockDrift = ([diff(cat(1,X(k).RecordStartTimeMATLAB));NaN]*86400 – cat(1,X(k).NumberSamples)./cat(1,X(k).SampleRate))./cat(1,X(k).NumberSamples);
I(i).OverlapBlockIndex = k(find(I(i).ClockDrift.*cat(1,X(k).NumberSamples).*cat(1,X(k).SampleRate) < -.5) + 1);
I(i).OverlapTime = cat(1,X(I(i).OverlapBlockIndex).RecordStartTimeMATLAB);
I(i).GapBlockIndex = k(find(I(i).ClockDrift.*cat(1,X(k).NumberSamples).*cat(1,X(k).SampleRate) > .5) + 1);
I(i).GapTime = cat(1,X(I(i).GapBlockIndex).RecordStartTimeMATLAB);
end
end
if nargout > 1
varargout{2} = I;
end
% — plots the data
if makeplot
figure
xlim = [min(cat(1,X.t)),max(cat(1,X.t))];
% test if all data records have the same length
rl = unique(cat(1,X.DataRecordSize));
if numel(rl) == 1
rl_text = sprintf(‘%d bytes’,rl);
else
rl_text = sprintf(‘%d-%d bytes’,min(rl),max(rl));
end
% test if all data records have the same sampling rate
sr = unique(cat(1,X.SampleRate));
if numel(sr) == 1
sr_text = sprintf(‘%g Hz’,sr);
else
sr_text = sprintf(‘%d # samp. rates’,numel(sr));
end
% test if all data records have the same encoding format
ef = unique(cellstr(cat(1,X.EncodingFormatName)));
if numel(ef) == 1
ef_text = sprintf(‘%s’,ef{:});
else
ef_text = sprintf(‘%d different encod. formats’,numel(ef));
end
if nc == 1
plot(cat(1,X.t),cat(1,X.d))
hold on
for i = 1:length(I.GapBlockIndex)
plot(I.GapTime(i),X(I.GapBlockIndex(i)).d(1),’*r’)
end
for i = 1:length(I.OverlapBlockIndex)
plot(I.OverlapTime(i),X(I.OverlapBlockIndex(i)).d(1),’og’)
end
hold off
set(gca,’XLim’,xlim)
datetick(‘x’,’keeplimits’)
grid on
xlabel(sprintf(‘Timen(%s to %s)’,datestr(xlim(1)),datestr(xlim(2))))
ylabel(‘Counts’)
title(sprintf(‘mini-SEED file "%s"n%s (%d rec. @ %s – %g samp. @ %s – %s)’, …
f,un{1},length(X),rl_text,numel(cat(1,X.d)),sr_text,ef_text),’Interpreter’,’none’)
else
% plot is done only for real data channels…
if nc > max_channels
warning(‘Plot has been limited to %d channels (over %d). See help to manage multiplexed file.’, …
max_channels,nc);
nc = max_channels;
end
for i = 1:nc
subplot(nc*2,1,i*2 + (-1:0))
k = I(i).XBlockIndex;
if ~any(strcmp(‘ASCII’,cellstr(cat(1,X(k).EncodingFormatName))))
plot(cat(1,X(k).t),cat(1,X(k).d))
hold on
for ii = 1:length(I(i).GapBlockIndex)
if ~isempty(X(I(i).GapBlockIndex(ii)).d)
plot(I(i).GapTime(ii),X(I(i).GapBlockIndex(ii)).d,’r’)
else
plot(repmat(I(i).GapTime(ii),1,2),ylim,’r’)
end
end
for ii = 1:length(I(i).OverlapBlockIndex)
if ~isempty(X(I(i).OverlapBlockIndex(ii)).d)
plot(I(i).OverlapTime(ii),X(I(i).OverlapBlockIndex(ii)).d,’g’)
else
plot(repmat(I(i).OverlapTime(ii),1,2),ylim,’g’)
end
end
hold off
end
set(gca,’XLim’,xlim,’FontSize’,8)
h = ylabel(un{i},’Interpreter’,’none’);
if nc > max_channel_label
set(gca,’YTick’,[])
set(h,’Rotation’,0,’HorizontalAlignment’,’right’,’FontSize’,8)
end
datetick(‘x’,’keeplimits’)
set(gca,’XTickLabel’,[])
grid on
if i == 1
title(sprintf(‘mini-SEED file "%s"n%d channels (%d rec. @ %s – %g data – %s – %s)’, …
f,length(un),length(X),rl_text,numel(cat(1,X(k).d)),sr_text,ef_text),’Interpreter’,’none’)
end
if i == nc
datetick(‘x’,’keeplimits’)
xlabel(sprintf(‘Timen(%s to %s)’,datestr(xlim(1)),datestr(xlim(2))))
end
end
v = version;
if str2double(v(1))>=7
linkaxes(findobj(gcf,’type’,’axes’),’x’)
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function D = read_data_record
% read_data_record uses global variables f, fid, offset, le, ef, wo, rl,
% and verbose. It reads a data record and returns a structure D.
global f fid offset le ef wo rl verbose notc force
fseek(fid,offset,’bof’);
% — read fixed section of Data Header (48 bytes)
D.SequenceNumber = fread(fid,6,’*char’)’;
D.DataQualityIndicator = fread(fid,1,’*char’);
D.ReservedByte = fread(fid,1,’*char’);
D.StationIdentifierCode = fread(fid,5,’*char’)’;
D.LocationIdentifier = fread(fid,2,’*char’)’;
D.ChannelIdentifier = fread(fid,3,’*char’)’;
D.NetworkCode = fread(fid,2,’*char’)’;
D.ChannelFullName = sprintf(‘%s:%s:%s:%s’,deblank(D.NetworkCode), …
deblank(D.StationIdentifierCode),deblank(D.LocationIdentifier), …
deblank(D.ChannelIdentifier));
% Start Time decoding
[D.RecordStartTime,swapflag] = readbtime;
D.RecordStartTimeISO = sprintf(‘%4d-%03d %02d:%02d:%07.4f’,D.RecordStartTime);
if swapflag
if le
machinefmt = ‘ieee-be’;
le = 0;
else
machinefmt = ‘ieee-le’;
le = 1;
end
position = ftell(fid);
fclose(fid);
fid = fopen(f,’rb’,machinefmt);
fseek(fid,position,’bof’);
if verbose > 0
warning(‘RDMSEED:DataIntegrity’, …
‘Sequence # %s: need to switch file encoding to %s…n’, …
D.SequenceNumber,machinefmt);
end
end
D.NumberSamples = fread(fid,1,’uint16′);
% Sample Rate decoding
SampleRateFactor = fread(fid,1,’int16′);
SampleRateMultiplier = fread(fid,1,’int16′);
if SampleRateFactor > 0
if SampleRateMultiplier >= 0
D.SampleRate = SampleRateFactor*SampleRateMultiplier;
else
D.SampleRate = -1*SampleRateFactor/SampleRateMultiplier;
end
else
if SampleRateMultiplier >= 0
D.SampleRate = -1*SampleRateMultiplier/SampleRateFactor;
else
D.SampleRate = 1/(SampleRateFactor*SampleRateMultiplier);
end
end
D.ActivityFlags = fread(fid,1,’uint8′);
D.IOFlags = fread(fid,1,’uint8′);
D.DataQualityFlags = fread(fid,1,’uint8′);
D.NumberBlockettesFollow = fread(fid,1,’uint8′);
D.TimeCorrection = fread(fid,1,’int32′); % Time correction in 0.0001 s
D.OffsetBeginData = fread(fid,1,’uint16′);
D.OffsetFirstBlockette = fread(fid,1,’uint16′);
% — read the blockettes
OffsetNextBlockette = D.OffsetFirstBlockette;
D.BLOCKETTES = [];
b2000 = 0; % Number of Blockette 2000
for i = 1:D.NumberBlockettesFollow
fseek(fid,offset + OffsetNextBlockette,’bof’);
BlocketteType = fread(fid,1,’uint16′);
switch BlocketteType
case 1000
% BLOCKETTE 1000 = Data Only SEED (8 bytes)
OffsetNextBlockette = fread(fid,1,’uint16′);
D.BLOCKETTES.B1000.EncodingFormat = fread(fid,1,’uint8′);
D.BLOCKETTES.B1000.WordOrder = fread(fid,1,’uint8′);
D.BLOCKETTES.B1000.DataRecordLength = fread(fid,1,’uint8′);
D.BLOCKETTES.B1000.Reserved = fread(fid,1,’uint8′);
case 1001
% BLOCKETTE 1001 = Data Extension (8 bytes)
OffsetNextBlockette = fread(fid,1,’uint16′);
D.BLOCKETTES.B1001.TimingQuality = fread(fid,1,’uint8′);
D.BLOCKETTES.B1001.Micro_sec = fread(fid,1,’int8′);
D.BLOCKETTES.B1001.Reserved = fread(fid,1,’uint8′);
D.BLOCKETTES.B1001.FrameCount = fread(fid,1,’uint8′);
case 100
% BLOCKETTE 100 = Sample Rate (12 bytes)
OffsetNextBlockette = fread(fid,1,’uint16′);
D.BLOCKETTES.B100.ActualSampleRate = fread(fid,1,’float32′);
D.BLOCKETTES.B100.Flags = fread(fid,1,’uint8′);
D.BLOCKETTES.B100.Reserved = fread(fid,1,’uint8′);
case 500
% BLOCKETTE 500 = Timing (200 bytes)
OffsetNextBlockette = fread(fid,1,’uint16′);
D.BLOCKETTES.B500.VCOCorrection = fread(fid,1,’float32′);
D.BLOCKETTES.B500.TimeOfException = readbtime;
D.BLOCKETTES.B500.MicroSec = fread(fid,1,’int8′);
D.BLOCKETTES.B500.ReceptionQuality = fread(fid,1,’uint8′);
D.BLOCKETTES.B500.ExceptionCount = fread(fid,1,’uint16′);
D.BLOCKETTES.B500.ExceptionType = fread(fid,16,’*char’)’;
D.BLOCKETTES.B500.ClockModel = fread(fid,32,’*char’)’;
D.BLOCKETTES.B500.ClockStatus = fread(fid,128,’*char’)’;
case 2000
% BLOCKETTE 2000 = Opaque Data (variable length)
b2000 = b2000 + 1;
OffsetNextBlockette = fread(fid,1,’uint16′);
BlocketteLength = fread(fid,1,’uint16′);
OffsetOpaqueData = fread(fid,1,’uint16′);
D.BLOCKETTES.B2000(b2000).RecordNumber = fread(fid,1,’uint32′);
D.BLOCKETTES.B2000(b2000).DataWordOrder = fread(fid,1,’uint8′);
D.BLOCKETTES.B2000(b2000).Flags = fread(fid,1,’uint8′);
NumberHeaderFields = fread(fid,1,’uint8′);
HeaderFields = splitfield(fread(fid,OffsetOpaqueData-15,’*char’)’,’~’);
D.BLOCKETTES.B2000(b2000).HeaderFields = HeaderFields(1:NumberHeaderFields);
% Opaque data are stored as a single char string, but must be
% decoded using appropriate format (e.g., Quanterra Q330)
D.BLOCKETTES.B2000(b2000).OpaqueData = fread(fid,BlocketteLength-OffsetOpaqueData,’*char’)’;
otherwise
OffsetNextBlockette = fread(fid,1,’uint16′);
if verbose > 0
warning(‘RDMSEED:UnknownBlockette’, …
‘Unknown Blockette number %d (%s)!n’, …
BlocketteType,D.ChannelFullName);
end
end
end
% — read the data stream
fseek(fid,offset + D.OffsetBeginData,’bof’);
if ~force && isfield(D.BLOCKETTES,’B1000′)
EncodingFormat = D.BLOCKETTES.B1000.EncodingFormat;
WordOrder = D.BLOCKETTES.B1000.WordOrder;
D.DataRecordSize = 2^D.BLOCKETTES.B1000.DataRecordLength;
else
EncodingFormat = ef;
WordOrder = wo;
D.DataRecordSize = rl;
end
uncoded = 0;
D.d = NaN;
D.t = NaN;
switch EncodingFormat
case 0
% — decoding format: ASCII text
D.EncodingFormatName = {‘ASCII’};
D.d = fread(fid,D.DataRecordSize – D.OffsetBeginData,’*char’)’;
case 1
% — decoding format: 16-bit integers
D.EncodingFormatName = {‘INT16’};
dd = fread(fid,ceil((D.DataRecordSize – D.OffsetBeginData)/2),’*int16′);
if xor(~WordOrder,le)
dd = swapbytes(dd);
end
D.d = dd(1:D.NumberSamples);
case 2
% — decoding format: 24-bit integers
D.EncodingFormatName = {‘INT24’};
dd = fread(fid,ceil((D.DataRecordSize – D.OffsetBeginData)/3),’bit24=>int32′);
if xor(~WordOrder,le)
dd = swapbytes(dd);
end
D.d = dd(1:D.NumberSamples);
case 3
% — decoding format: 32-bit integers
D.EncodingFormatName = {‘INT32’};
dd = fread(fid,ceil((D.DataRecordSize – D.OffsetBeginData)/4),’*int32′);
if xor(~WordOrder,le)
dd = swapbytes(dd);
end
D.d = dd(1:D.NumberSamples);
case 4
% — decoding format: IEEE floating point
D.EncodingFormatName = {‘FLOAT32’};
dd = fread(fid,ceil((D.DataRecordSize – D.OffsetBeginData)/4),’*float’);
if xor(~WordOrder,le)
dd = swapbytes(dd);
end
D.d = dd(1:D.NumberSamples);
case 5
% — decoding format: IEEE double precision floating point
D.EncodingFormatName = {‘FLOAT64’};
dd = fread(fid,ceil((D.DataRecordSize – D.OffsetBeginData)/8),’*double’);
if xor(~WordOrder,le)
dd = swapbytes(dd);
end
D.d = dd(1:D.NumberSamples);
case {10,11,19}
% — decoding formats: STEIM-1 and STEIM-2 compression
% (c) Joseph M. Steim, Quanterra Inc., 1994
steim = find(EncodingFormat==[10,11,19]);
D.EncodingFormatName = {sprintf(‘STEIM%d’,steim)};
% Steim compression decoding strategy optimized for Matlab
% — by F. Beauducel, October 2010 —
%
% 1. loads all data into a single 16xM uint32 array
% 2. gets all nibbles from the first row splitted into 2-bit values
% 3. for each possible nibble value, selects (find) and decodes
% (bitsplit) all the corresponding words, and stores results
% in a 4xN (STEIM1) or 7xN (STEIM2) array previously filled with
% NaN’s. For STEIM2 with nibbles 2 or 3, decodes also dnib values
% (first 2-bit of the word)
% 5. reduces this array with non-NaN values only
% 6. integrates with cumsum
%
% This method is about 30 times faster than a ‘C-like’ loops coding…
frame32 = fread(fid,[16,(D.DataRecordSize – D.OffsetBeginData)/64],’*uint32′);
if xor(~WordOrder,le)
frame32 = swapbytes(frame32);
end
% specific processes for STEIM-3
if steim == 3
% first bit = 1 means second differences
SecondDiff = bitshift(frame32(1,:),-31);
% checks for "squeezed flag"… and replaces frame32(1,:)
squeezed = bitand(bitshift(frame32(1,:),-24),127);
k = find(bitget(squeezed,7));
if ~isempty(k)
moredata24 = bitand(frame32(1,k),16777215);
k = find(squeezed == 80); % upper nibble 8-bit = 0x50
if ~isempty(k)
frame32(1,k) = hex2dec(‘15555555’);
end
k = find(squeezed == 96); % upper nibble 8-bit = 0x60
if ~isempty(k)
frame32(1,k) = hex2dec(‘2aaaaaaa’);
end
k = find(squeezed == 112); % upper nibble 8-bit = 0x70
if ~isempty(k)
frame32(1,k) = hex2dec(‘3fffffff’);
end
end
end
% nibbles is an array of the same size as frame32…
nibbles = bitand(bitshift(repmat(frame32(1,:),16,1),repmat(-30:2:0,size(frame32,2),1)’),3);
x0 = bitsign(frame32(2,1),32); % forward integration constant
xn = bitsign(frame32(3,1),32); % reverse integration constant
switch steim
case 1
% STEIM-1: 3 cases following the nibbles
ddd = NaN*ones(4,numel(frame32)); % initiates array with NaN
k = find(nibbles == 1); % nibble = 1 : four 8-bit differences
if ~isempty(k)
ddd(1:4,k) = bitsplit(frame32(k),32,8);
end
k = find(nibbles == 2); % nibble = 2 : two 16-bit differences
if ~isempty(k)
ddd(1:2,k) = bitsplit(frame32(k),32,16);
end
k = find(nibbles == 3); % nibble = 3 : one 32-bit difference
if ~isempty(k)
ddd(1,k) = bitsign(frame32(k),32);
end
case 2
% STEIM-2: 7 cases following the nibbles and dnib
ddd = NaN*ones(7,numel(frame32)); % initiates array with NaN
k = find(nibbles == 1); % nibble = 1 : four 8-bit differences
if ~isempty(k)
ddd(1:4,k) = bitsplit(frame32(k),32,8);
end
k = find(nibbles == 2); % nibble = 2 : must look in dnib
if ~isempty(k)
dnib = bitshift(frame32(k),-30);
kk = k(dnib == 1); % dnib = 1 : one 30-bit difference
if ~isempty(kk)
ddd(1,kk) = bitsign(frame32(kk),30);
end
kk = k(dnib == 2); % dnib = 2 : two 15-bit differences
if ~isempty(kk)
ddd(1:2,kk) = bitsplit(frame32(kk),30,15);
end
kk = k(dnib == 3); % dnib = 3 : three 10-bit differences
if ~isempty(kk)
ddd(1:3,kk) = bitsplit(frame32(kk),30,10);
end
end
k = find(nibbles == 3); % nibble = 3 : must look in dnib
if ~isempty(k)
dnib = bitshift(frame32(k),-30);
kk = k(dnib == 0); % dnib = 0 : five 6-bit difference
if ~isempty(kk)
ddd(1:5,kk) = bitsplit(frame32(kk),30,6);
end
kk = k(dnib == 1); % dnib = 1 : six 5-bit differences
if ~isempty(kk)
ddd(1:6,kk) = bitsplit(frame32(kk),30,5);
end
kk = k(dnib == 2); % dnib = 2 : seven 4-bit differences (28 bits!)
if ~isempty(kk)
ddd(1:7,kk) = bitsplit(frame32(kk),28,4);
end
end
case 3 % *** STEIM-3 DECODING IS ALPHA AND UNTESTED ***
% STEIM-3: 7 cases following the nibbles
ddd = NaN*ones(9,numel(frame32)); % initiates array with NaN
k = find(nibbles == 0); % nibble = 0 : two 16-bit differences
if ~isempty(k)
ddd(1:2,k) = bitsplit(frame32(k),32,16);
end
k = find(nibbles == 1); % nibble = 1 : four 8-bit differences
if ~isempty(k)
ddd(1:4,k) = bitsplit(frame32(k),32,8);
end
k = find(nibbles == 2); % nibble = 2 : must look even dnib
if ~isempty(k)
dnib2 = bitshift(frame32(k(2:2:end)),-30);
w60 = bitand(frame32(k(2:2:end)),1073741823) …
+ bitshift(bitand(frame32(k(1:2:end)),1073741823),30); % concatenates two 30-bit words
kk = find(dnib2 == 0); % dnib = 0: five 12-bit differences (60 bits)
if ~isempty(kk)
ddd(1:5,k(2*kk)) = bitsplit(w60,60,12);
end
kk = find(dnib2 == 1); % dnib = 1: three 20-bit differences (60 bits)
if ~isempty(kk)
ddd(1:3,k(2*kk)) = bitsplit(w60,60,20);
end
end
k = find(nibbles == 3); % nibble = 3 : must look 3rd bit
if ~isempty(k)
dnib = bitshift(frame32(k),-27);
kk = k(dnib == 24); % dnib = 11000 : nine 3-bit differences (27 bits)
if ~isempty(kk)
ddd(1:9,kk) = bitsplit(frame32(kk),27,3);
end
kk = k(dnib == 25); % dnib = 11001 : Not A Difference
if ~isempty(kk)
ddd(1,kk) = bitsign(frame32(kk),27);
end
kk = k(dnib > 27); % dnib = 111.. : 29-bit sample (29 bits)
if ~isempty(kk)
ddd(1,kk) = bitsign(frame32(kk),29);
end
end
end
% Little-endian coding: needs to swap bytes
if ~WordOrder
ddd = flipud(ddd);
end
dd = ddd(~isnan(ddd)); % reduces initial array ddd: dd is non-NaN values of ddd
% controls the number of samples
if numel(dd) ~= D.NumberSamples
if verbose > 1
warning(‘RDMSEED:DataIntegrity’,’Problem in %s sequence # %s [%d-%03d %02d:%02d:%07.4f]: number of samples in header (%d) does not equal data (%d).n’, …
D.EncodingFormatName{:},D.SequenceNumber,D.RecordStartTimeISO,D.NumberSamples,numel(dd));
end
if numel(dd) < D.NumberSamples
D.NumberSamples = numel(dd);
end
end
% rebuilds the data vector by integrating the differences
D.d = cumsum([x0;dd(2:D.NumberSamples)]);
% controls data integrity…
if D.d(end) ~= xn
warning(‘RDMSEED:DataIntegrity’,’Problem in %s sequence # %s [%s]: data integrity check failed, last_data=%d, Xn=%d.n’, …
D.EncodingFormatName{:},D.SequenceNumber,D.RecordStartTimeISO,D.d(end),xn);
end
if D.NumberSamples == 0
D.d = nan(0,1);
end
% for debug purpose…
if verbose > 2
D.dd = dd;
D.nibbles = nibbles;
D.x0 = x0;
D.xn = xn;
end
case 12
% — decoding format: GEOSCOPE multiplexed 24-bit integer
D.EncodingFormatName = {‘GEOSCOPE24’};
dd = fread(fid,(D.DataRecordSize – D.OffsetBeginData)/3,’bit24=>double’);
if xor(~WordOrder,le)
dd = swapbytes(dd);
end
D.d = dd(1:D.NumberSamples);
case {13,14}
% — decoding format: GEOSCOPE multiplexed 16/3 and 16/4 bit gain ranged
% (13): 16/3-bit (bit 15 is unused)
% (14): 16/4-bit
% bits 15-12 = 3 or 4-bit gain exponent (positive)
% bits 11-0 = 12-bit mantissa (positive)
% => data = (mantissa – 2048) / 2^gain
geoscope = 7 + 8*(EncodingFormat==14); % mask for gain exponent
D.EncodingFormatName = {sprintf(‘GEOSCOPE16-%d’,EncodingFormat-10)};
dd = fread(fid,(D.DataRecordSize – D.OffsetBeginData)/2,’*uint16′);
if xor(~WordOrder,le)
dd = swapbytes(dd);
end
dd = (double(bitand(dd,2^12-1))-2^11)./2.^double(bitand(bitshift(dd,-12),geoscope));
D.d = dd(1:D.NumberSamples);
case 15
% — decoding format: US National Network compression
D.EncodingFormatName = {‘USNN’};
uncoded = 1;
case 16
% — decoding format: CDSN 16-bit gain ranged
D.EncodingFormatName = {‘CDSN’};
uncoded = 1;
case 17
% — decoding format: Graefenberg 16-bit gain ranged
D.EncodingFormatName = {‘GRAEFENBERG’};
uncoded = 1;
case 18
% — decoding format: IPG – Strasbourg 16-bit gain ranged
D.EncodingFormatName = {‘IPGS’};
uncoded = 1;
case 30
% — decoding format: SRO format
D.EncodingFormatName = {‘SRO’};
uncoded = 1;
case 31
% — decoding format: HGLP format
D.EncodingFormatName = {‘HGLP’};
uncoded = 1;
case 32
% — decoding format: DWWSSN gain ranged format
D.EncodingFormatName = {‘DWWSSN’};
uncoded = 1;
case 33
% — decoding format: RSTN 16-bit gain ranged
D.EncodingFormatName = {‘RSTN’};
uncoded = 1;
otherwise
D.EncodingFormatName = {sprintf(‘** Unknown (%d) **’,EncodingFormat)};
uncoded = 1;
end
if uncoded
error(‘Sorry, the encoding format "%s" is not yet implemented.’,D.EncodingFormatName);
end
% Applies time correction (if needed)
D.RecordStartTimeMATLAB = datenum(double([D.RecordStartTime(1),0,D.RecordStartTime(2:5)])) …
+ (~notc & bitand(D.ActivityFlags,2) == 0)*D.TimeCorrection/1e4/86400;
tv = datevec(D.RecordStartTimeMATLAB);
doy = datenum(tv(1:3)) – datenum(tv(1),1,0);
D.RecordStartTime = [tv(1),doy,tv(4:5),round(tv(6)*1e4)/1e4];
D.RecordStartTimeISO = sprintf(‘%4d-%03d %02d:%02d:%07.4f’,D.RecordStartTime);
D.t = D.RecordStartTimeMATLAB;
% makes the time vector and applies time correction (if needed)
if EncodingFormat > 0
D.t = D.t + (0:(D.NumberSamples-1))’/(D.SampleRate*86400);
end
offset = ftell(fid);
fread(fid,1,’char’); % this is to force EOF=1 on last record.
if feof(fid)
offset = -1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function c = splitfield(s,d)
% splitfield(S) splits string S of D-character separated field names
C = textscan(s,’%s’,’Delimiter’,d);
c = C{1};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [d,swapflag] = readbtime
% readbtime reads BTIME structure from current opened file and returns
% D = [YEAR,DAY,HOUR,MINUTE,SECONDS]
global fid forcebe
Year = fread(fid,1,’*uint16′);
DayOfYear = fread(fid,1,’*uint16′);
Hours = fread(fid,1,’uint8′);
Minutes = fread(fid,1,’uint8′);
Seconds = fread(fid,1,’uint8′);
fseek(fid,1,0); % skip 1 byte (unused)
Seconds0001 = fread(fid,1,’*uint16′);
% Automatic detection of little/big-endian encoding
% — by F. Beauducel, March 2014 —
%
% If the 2-byte day is >= 512, the file is not opened in the correct
% endianness. If the day is 1 or 256, there is a possible byte-swap and we
% need to check also the year; but we need to consider what is a valid year:
% – years from 1801 to 2047 are OK (swapbytes >= 2312)
% – years from 2048 to 2055 are OK (swapbytes <= 1800)
% – year 2056 is ambiguous (swapbytes = 2056)
% – years from 2057 to 2311 are OK (swapbytes >= 2312)
% – year 1799 is ambiguous (swapbytes = 1799)
% – year 1800 is suspicious (swapbytes = 2055)
%
% Thus, the only cases for which we are ‘sure’ there is a byte-swap, are:
% – day >= 512
% – (day == 1 or day == 256) and (year < 1799 or year > 2311)
%
% Note: in IRIS libmseed, the test is only year>2050 or year<1920.
if ~forcebe && (DayOfYear >= 512 || (ismember(DayOfYear,[1,256]) && (Year > 2311 || Year < 1799)))
swapflag = 1;
Year = swapbytes(Year);
DayOfYear = swapbytes(DayOfYear);
Seconds0001 = swapbytes(Seconds0001);
else
swapflag = 0;
end
d = [double(Year),double(DayOfYear),Hours,Minutes,Seconds + double(Seconds0001)/1e4];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function d = bitsplit(x,b,n)
% bitsplit(X,B,N) splits the B-bit number X into signed N-bit array
% X must be unsigned integer class
% N ranges from 1 to B
% B is a multiple of N
sign = repmat((b:-n:n)’,1,size(x,1));
x = repmat(x’,b/n,1);
d = double(bitand(bitshift(x,flipud(sign-b)),2^n-1)) …
– double(bitget(x,sign))*2^n;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function d = bitsign(x,n)
% bitsign(X,N) returns signed double value from unsigned N-bit number X.
% This is equivalent to bitsplit(X,N,N), but the formula is simplified so
% it is much more efficient
d = double(bitand(x,2^n-1)) – double(bitget(x,n)).*2^n;How can i change the subplots sequence in rdmseed.m, while plotting a miniseed file with e.g. 5 stations?
function varargout = rdmseed(varargin)
if nargin > 6
error(‘Too many input arguments.’)
end
% global variables shared with sub-functions
global f fid offset le ef wo rl forcebe verbose notc force
% default input arguments
makeplot = 0; % make plot flag
verbose = 0; % verbose flag/level
forcebe = 0; % force big-endian
ef = 10; % encoding format default
wo = 1; % word order default
rl = 2^12; % record length default
force = 0; % force input argument over blockette 1000 (UNDOCUMENTED)
notc = 0; % force no time correction (over ActivityFlags)
nullhead = 0; % allow null bytes before header
if nargin < 1
[filename,pathname] = uigetfile(‘*’,’Please select a miniSEED file…’);
f = fullfile(pathname,filename);
else
f = varargin{1};
end
if ~ischar(f) || ~exist(f,’file’)
error(‘File %s does not exist.’,f);
end
if nargin > 1
verbose = any(strcmpi(varargin,’v’)) + 2*any(strcmpi(varargin,’vv’)) …
+ 3*any(strcmpi(varargin,’vvv’));
makeplot = any(strcmpi(varargin,’plot’));
forcebe = any(strcmpi(varargin,’be’));
notc = any(strcmpi(varargin,’notc’));
force = any(strcmpi(varargin,’force’));
nullhead = any(strcmpi(varargin,’nullhead’));
end
nargs = (makeplot>0) + (verbose>0) + (forcebe>0) + (notc>0) + (force>0) …
+ (nullhead>0);
if nargin > (1 + nargs)
ef = varargin{2};
if ~isnumeric(ef) || ~any(ef==[0:5,10:19,30:33])
error(‘Argument ENCODINGFORMAT must be a valid FDSN code value.’);
end
end
if nargin > (2 + nargs)
wo = varargin{3};
if ~isnumeric(wo) || (wo ~= 0 && wo ~= 1)
error(‘Argument WORDORDER must be 0 or 1.’);
end
end
if nargin > (3 + nargs)
rl = varargin{4};
if ~isnumeric(rl) || rl < 256 || rem(log(rl)/log(2),1) ~= 0
error(‘Argument RECORDLENGTH must be a power of 2 and greater or equal to 256.’);
end
end
if nargout == 0
makeplot = 1;
end
% sensible limits for multiplexed files
max_channels = 20; % absolute max number of channels to plot
max_channel_label = 6; % max. number of channels for y-labels
% file is opened in Big-Endian encoding (this is encouraged by SEED)
fid = fopen(f,’rb’,’ieee-be’);
le = 0;
offset = 0;
% — tests if the header is mini-SEED
% the 7th character must be one of the "data header/quality indicator", usually ‘D’
header = fread(fid,20,’*char’);
if ~ismember(header(7),’DRMQ’)
if ismember(header(7),’VAST’)
error(‘File seems to be a SEED Volume. Cannot read it.’);
else
if header(1)==0
if nullhead
if verbose
fprintf(‘Null header option: bypassing…’);
end
c = 0;
fseek(fid,0,’bof’);
while c==0
c = fread(fid,1,’*char’);
offset = offset + 1;
end
if verbose
fprintf(‘ %d null bytes.n’,offset);
end
header = fread(fid,6,’*char’);
if ~ismember(header(6),’DRMQ’)
error(‘File is not in mini-SEED format. Cannot read it.’);
else
offset = offset – 1;
end
else
error(‘File starts with null bytes… if you believe it is still a miniseed file, try the ”nullhead” option.’);
end
else
error(‘File is not in mini-SEED format. Cannot read it.’);
end
end
end
i = 1;
while offset >= 0
X(i) = read_data_record;
i = i + 1;
end
fclose(fid);
if nargout > 0
varargout{1} = X;
end
% — analyses data
if makeplot || nargout > 1
% test if the file is multiplexed or a single channel
un = unique(cellstr(char(X.ChannelFullName)));
nc = numel(un);
for i = 1:nc
k = find(strcmp(cellstr(char(X.ChannelFullName)),un{i}));
I(i).ChannelFullName = X(k(1)).ChannelFullName;
I(i).XBlockIndex = k;
I(i).ClockDrift = ([diff(cat(1,X(k).RecordStartTimeMATLAB));NaN]*86400 – cat(1,X(k).NumberSamples)./cat(1,X(k).SampleRate))./cat(1,X(k).NumberSamples);
I(i).OverlapBlockIndex = k(find(I(i).ClockDrift.*cat(1,X(k).NumberSamples).*cat(1,X(k).SampleRate) < -.5) + 1);
I(i).OverlapTime = cat(1,X(I(i).OverlapBlockIndex).RecordStartTimeMATLAB);
I(i).GapBlockIndex = k(find(I(i).ClockDrift.*cat(1,X(k).NumberSamples).*cat(1,X(k).SampleRate) > .5) + 1);
I(i).GapTime = cat(1,X(I(i).GapBlockIndex).RecordStartTimeMATLAB);
end
end
if nargout > 1
varargout{2} = I;
end
% — plots the data
if makeplot
figure
xlim = [min(cat(1,X.t)),max(cat(1,X.t))];
% test if all data records have the same length
rl = unique(cat(1,X.DataRecordSize));
if numel(rl) == 1
rl_text = sprintf(‘%d bytes’,rl);
else
rl_text = sprintf(‘%d-%d bytes’,min(rl),max(rl));
end
% test if all data records have the same sampling rate
sr = unique(cat(1,X.SampleRate));
if numel(sr) == 1
sr_text = sprintf(‘%g Hz’,sr);
else
sr_text = sprintf(‘%d # samp. rates’,numel(sr));
end
% test if all data records have the same encoding format
ef = unique(cellstr(cat(1,X.EncodingFormatName)));
if numel(ef) == 1
ef_text = sprintf(‘%s’,ef{:});
else
ef_text = sprintf(‘%d different encod. formats’,numel(ef));
end
if nc == 1
plot(cat(1,X.t),cat(1,X.d))
hold on
for i = 1:length(I.GapBlockIndex)
plot(I.GapTime(i),X(I.GapBlockIndex(i)).d(1),’*r’)
end
for i = 1:length(I.OverlapBlockIndex)
plot(I.OverlapTime(i),X(I.OverlapBlockIndex(i)).d(1),’og’)
end
hold off
set(gca,’XLim’,xlim)
datetick(‘x’,’keeplimits’)
grid on
xlabel(sprintf(‘Timen(%s to %s)’,datestr(xlim(1)),datestr(xlim(2))))
ylabel(‘Counts’)
title(sprintf(‘mini-SEED file "%s"n%s (%d rec. @ %s – %g samp. @ %s – %s)’, …
f,un{1},length(X),rl_text,numel(cat(1,X.d)),sr_text,ef_text),’Interpreter’,’none’)
else
% plot is done only for real data channels…
if nc > max_channels
warning(‘Plot has been limited to %d channels (over %d). See help to manage multiplexed file.’, …
max_channels,nc);
nc = max_channels;
end
for i = 1:nc
subplot(nc*2,1,i*2 + (-1:0))
k = I(i).XBlockIndex;
if ~any(strcmp(‘ASCII’,cellstr(cat(1,X(k).EncodingFormatName))))
plot(cat(1,X(k).t),cat(1,X(k).d))
hold on
for ii = 1:length(I(i).GapBlockIndex)
if ~isempty(X(I(i).GapBlockIndex(ii)).d)
plot(I(i).GapTime(ii),X(I(i).GapBlockIndex(ii)).d,’r’)
else
plot(repmat(I(i).GapTime(ii),1,2),ylim,’r’)
end
end
for ii = 1:length(I(i).OverlapBlockIndex)
if ~isempty(X(I(i).OverlapBlockIndex(ii)).d)
plot(I(i).OverlapTime(ii),X(I(i).OverlapBlockIndex(ii)).d,’g’)
else
plot(repmat(I(i).OverlapTime(ii),1,2),ylim,’g’)
end
end
hold off
end
set(gca,’XLim’,xlim,’FontSize’,8)
h = ylabel(un{i},’Interpreter’,’none’);
if nc > max_channel_label
set(gca,’YTick’,[])
set(h,’Rotation’,0,’HorizontalAlignment’,’right’,’FontSize’,8)
end
datetick(‘x’,’keeplimits’)
set(gca,’XTickLabel’,[])
grid on
if i == 1
title(sprintf(‘mini-SEED file "%s"n%d channels (%d rec. @ %s – %g data – %s – %s)’, …
f,length(un),length(X),rl_text,numel(cat(1,X(k).d)),sr_text,ef_text),’Interpreter’,’none’)
end
if i == nc
datetick(‘x’,’keeplimits’)
xlabel(sprintf(‘Timen(%s to %s)’,datestr(xlim(1)),datestr(xlim(2))))
end
end
v = version;
if str2double(v(1))>=7
linkaxes(findobj(gcf,’type’,’axes’),’x’)
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function D = read_data_record
% read_data_record uses global variables f, fid, offset, le, ef, wo, rl,
% and verbose. It reads a data record and returns a structure D.
global f fid offset le ef wo rl verbose notc force
fseek(fid,offset,’bof’);
% — read fixed section of Data Header (48 bytes)
D.SequenceNumber = fread(fid,6,’*char’)’;
D.DataQualityIndicator = fread(fid,1,’*char’);
D.ReservedByte = fread(fid,1,’*char’);
D.StationIdentifierCode = fread(fid,5,’*char’)’;
D.LocationIdentifier = fread(fid,2,’*char’)’;
D.ChannelIdentifier = fread(fid,3,’*char’)’;
D.NetworkCode = fread(fid,2,’*char’)’;
D.ChannelFullName = sprintf(‘%s:%s:%s:%s’,deblank(D.NetworkCode), …
deblank(D.StationIdentifierCode),deblank(D.LocationIdentifier), …
deblank(D.ChannelIdentifier));
% Start Time decoding
[D.RecordStartTime,swapflag] = readbtime;
D.RecordStartTimeISO = sprintf(‘%4d-%03d %02d:%02d:%07.4f’,D.RecordStartTime);
if swapflag
if le
machinefmt = ‘ieee-be’;
le = 0;
else
machinefmt = ‘ieee-le’;
le = 1;
end
position = ftell(fid);
fclose(fid);
fid = fopen(f,’rb’,machinefmt);
fseek(fid,position,’bof’);
if verbose > 0
warning(‘RDMSEED:DataIntegrity’, …
‘Sequence # %s: need to switch file encoding to %s…n’, …
D.SequenceNumber,machinefmt);
end
end
D.NumberSamples = fread(fid,1,’uint16′);
% Sample Rate decoding
SampleRateFactor = fread(fid,1,’int16′);
SampleRateMultiplier = fread(fid,1,’int16′);
if SampleRateFactor > 0
if SampleRateMultiplier >= 0
D.SampleRate = SampleRateFactor*SampleRateMultiplier;
else
D.SampleRate = -1*SampleRateFactor/SampleRateMultiplier;
end
else
if SampleRateMultiplier >= 0
D.SampleRate = -1*SampleRateMultiplier/SampleRateFactor;
else
D.SampleRate = 1/(SampleRateFactor*SampleRateMultiplier);
end
end
D.ActivityFlags = fread(fid,1,’uint8′);
D.IOFlags = fread(fid,1,’uint8′);
D.DataQualityFlags = fread(fid,1,’uint8′);
D.NumberBlockettesFollow = fread(fid,1,’uint8′);
D.TimeCorrection = fread(fid,1,’int32′); % Time correction in 0.0001 s
D.OffsetBeginData = fread(fid,1,’uint16′);
D.OffsetFirstBlockette = fread(fid,1,’uint16′);
% — read the blockettes
OffsetNextBlockette = D.OffsetFirstBlockette;
D.BLOCKETTES = [];
b2000 = 0; % Number of Blockette 2000
for i = 1:D.NumberBlockettesFollow
fseek(fid,offset + OffsetNextBlockette,’bof’);
BlocketteType = fread(fid,1,’uint16′);
switch BlocketteType
case 1000
% BLOCKETTE 1000 = Data Only SEED (8 bytes)
OffsetNextBlockette = fread(fid,1,’uint16′);
D.BLOCKETTES.B1000.EncodingFormat = fread(fid,1,’uint8′);
D.BLOCKETTES.B1000.WordOrder = fread(fid,1,’uint8′);
D.BLOCKETTES.B1000.DataRecordLength = fread(fid,1,’uint8′);
D.BLOCKETTES.B1000.Reserved = fread(fid,1,’uint8′);
case 1001
% BLOCKETTE 1001 = Data Extension (8 bytes)
OffsetNextBlockette = fread(fid,1,’uint16′);
D.BLOCKETTES.B1001.TimingQuality = fread(fid,1,’uint8′);
D.BLOCKETTES.B1001.Micro_sec = fread(fid,1,’int8′);
D.BLOCKETTES.B1001.Reserved = fread(fid,1,’uint8′);
D.BLOCKETTES.B1001.FrameCount = fread(fid,1,’uint8′);
case 100
% BLOCKETTE 100 = Sample Rate (12 bytes)
OffsetNextBlockette = fread(fid,1,’uint16′);
D.BLOCKETTES.B100.ActualSampleRate = fread(fid,1,’float32′);
D.BLOCKETTES.B100.Flags = fread(fid,1,’uint8′);
D.BLOCKETTES.B100.Reserved = fread(fid,1,’uint8′);
case 500
% BLOCKETTE 500 = Timing (200 bytes)
OffsetNextBlockette = fread(fid,1,’uint16′);
D.BLOCKETTES.B500.VCOCorrection = fread(fid,1,’float32′);
D.BLOCKETTES.B500.TimeOfException = readbtime;
D.BLOCKETTES.B500.MicroSec = fread(fid,1,’int8′);
D.BLOCKETTES.B500.ReceptionQuality = fread(fid,1,’uint8′);
D.BLOCKETTES.B500.ExceptionCount = fread(fid,1,’uint16′);
D.BLOCKETTES.B500.ExceptionType = fread(fid,16,’*char’)’;
D.BLOCKETTES.B500.ClockModel = fread(fid,32,’*char’)’;
D.BLOCKETTES.B500.ClockStatus = fread(fid,128,’*char’)’;
case 2000
% BLOCKETTE 2000 = Opaque Data (variable length)
b2000 = b2000 + 1;
OffsetNextBlockette = fread(fid,1,’uint16′);
BlocketteLength = fread(fid,1,’uint16′);
OffsetOpaqueData = fread(fid,1,’uint16′);
D.BLOCKETTES.B2000(b2000).RecordNumber = fread(fid,1,’uint32′);
D.BLOCKETTES.B2000(b2000).DataWordOrder = fread(fid,1,’uint8′);
D.BLOCKETTES.B2000(b2000).Flags = fread(fid,1,’uint8′);
NumberHeaderFields = fread(fid,1,’uint8′);
HeaderFields = splitfield(fread(fid,OffsetOpaqueData-15,’*char’)’,’~’);
D.BLOCKETTES.B2000(b2000).HeaderFields = HeaderFields(1:NumberHeaderFields);
% Opaque data are stored as a single char string, but must be
% decoded using appropriate format (e.g., Quanterra Q330)
D.BLOCKETTES.B2000(b2000).OpaqueData = fread(fid,BlocketteLength-OffsetOpaqueData,’*char’)’;
otherwise
OffsetNextBlockette = fread(fid,1,’uint16′);
if verbose > 0
warning(‘RDMSEED:UnknownBlockette’, …
‘Unknown Blockette number %d (%s)!n’, …
BlocketteType,D.ChannelFullName);
end
end
end
% — read the data stream
fseek(fid,offset + D.OffsetBeginData,’bof’);
if ~force && isfield(D.BLOCKETTES,’B1000′)
EncodingFormat = D.BLOCKETTES.B1000.EncodingFormat;
WordOrder = D.BLOCKETTES.B1000.WordOrder;
D.DataRecordSize = 2^D.BLOCKETTES.B1000.DataRecordLength;
else
EncodingFormat = ef;
WordOrder = wo;
D.DataRecordSize = rl;
end
uncoded = 0;
D.d = NaN;
D.t = NaN;
switch EncodingFormat
case 0
% — decoding format: ASCII text
D.EncodingFormatName = {‘ASCII’};
D.d = fread(fid,D.DataRecordSize – D.OffsetBeginData,’*char’)’;
case 1
% — decoding format: 16-bit integers
D.EncodingFormatName = {‘INT16’};
dd = fread(fid,ceil((D.DataRecordSize – D.OffsetBeginData)/2),’*int16′);
if xor(~WordOrder,le)
dd = swapbytes(dd);
end
D.d = dd(1:D.NumberSamples);
case 2
% — decoding format: 24-bit integers
D.EncodingFormatName = {‘INT24’};
dd = fread(fid,ceil((D.DataRecordSize – D.OffsetBeginData)/3),’bit24=>int32′);
if xor(~WordOrder,le)
dd = swapbytes(dd);
end
D.d = dd(1:D.NumberSamples);
case 3
% — decoding format: 32-bit integers
D.EncodingFormatName = {‘INT32’};
dd = fread(fid,ceil((D.DataRecordSize – D.OffsetBeginData)/4),’*int32′);
if xor(~WordOrder,le)
dd = swapbytes(dd);
end
D.d = dd(1:D.NumberSamples);
case 4
% — decoding format: IEEE floating point
D.EncodingFormatName = {‘FLOAT32’};
dd = fread(fid,ceil((D.DataRecordSize – D.OffsetBeginData)/4),’*float’);
if xor(~WordOrder,le)
dd = swapbytes(dd);
end
D.d = dd(1:D.NumberSamples);
case 5
% — decoding format: IEEE double precision floating point
D.EncodingFormatName = {‘FLOAT64’};
dd = fread(fid,ceil((D.DataRecordSize – D.OffsetBeginData)/8),’*double’);
if xor(~WordOrder,le)
dd = swapbytes(dd);
end
D.d = dd(1:D.NumberSamples);
case {10,11,19}
% — decoding formats: STEIM-1 and STEIM-2 compression
% (c) Joseph M. Steim, Quanterra Inc., 1994
steim = find(EncodingFormat==[10,11,19]);
D.EncodingFormatName = {sprintf(‘STEIM%d’,steim)};
% Steim compression decoding strategy optimized for Matlab
% — by F. Beauducel, October 2010 —
%
% 1. loads all data into a single 16xM uint32 array
% 2. gets all nibbles from the first row splitted into 2-bit values
% 3. for each possible nibble value, selects (find) and decodes
% (bitsplit) all the corresponding words, and stores results
% in a 4xN (STEIM1) or 7xN (STEIM2) array previously filled with
% NaN’s. For STEIM2 with nibbles 2 or 3, decodes also dnib values
% (first 2-bit of the word)
% 5. reduces this array with non-NaN values only
% 6. integrates with cumsum
%
% This method is about 30 times faster than a ‘C-like’ loops coding…
frame32 = fread(fid,[16,(D.DataRecordSize – D.OffsetBeginData)/64],’*uint32′);
if xor(~WordOrder,le)
frame32 = swapbytes(frame32);
end
% specific processes for STEIM-3
if steim == 3
% first bit = 1 means second differences
SecondDiff = bitshift(frame32(1,:),-31);
% checks for "squeezed flag"… and replaces frame32(1,:)
squeezed = bitand(bitshift(frame32(1,:),-24),127);
k = find(bitget(squeezed,7));
if ~isempty(k)
moredata24 = bitand(frame32(1,k),16777215);
k = find(squeezed == 80); % upper nibble 8-bit = 0x50
if ~isempty(k)
frame32(1,k) = hex2dec(‘15555555’);
end
k = find(squeezed == 96); % upper nibble 8-bit = 0x60
if ~isempty(k)
frame32(1,k) = hex2dec(‘2aaaaaaa’);
end
k = find(squeezed == 112); % upper nibble 8-bit = 0x70
if ~isempty(k)
frame32(1,k) = hex2dec(‘3fffffff’);
end
end
end
% nibbles is an array of the same size as frame32…
nibbles = bitand(bitshift(repmat(frame32(1,:),16,1),repmat(-30:2:0,size(frame32,2),1)’),3);
x0 = bitsign(frame32(2,1),32); % forward integration constant
xn = bitsign(frame32(3,1),32); % reverse integration constant
switch steim
case 1
% STEIM-1: 3 cases following the nibbles
ddd = NaN*ones(4,numel(frame32)); % initiates array with NaN
k = find(nibbles == 1); % nibble = 1 : four 8-bit differences
if ~isempty(k)
ddd(1:4,k) = bitsplit(frame32(k),32,8);
end
k = find(nibbles == 2); % nibble = 2 : two 16-bit differences
if ~isempty(k)
ddd(1:2,k) = bitsplit(frame32(k),32,16);
end
k = find(nibbles == 3); % nibble = 3 : one 32-bit difference
if ~isempty(k)
ddd(1,k) = bitsign(frame32(k),32);
end
case 2
% STEIM-2: 7 cases following the nibbles and dnib
ddd = NaN*ones(7,numel(frame32)); % initiates array with NaN
k = find(nibbles == 1); % nibble = 1 : four 8-bit differences
if ~isempty(k)
ddd(1:4,k) = bitsplit(frame32(k),32,8);
end
k = find(nibbles == 2); % nibble = 2 : must look in dnib
if ~isempty(k)
dnib = bitshift(frame32(k),-30);
kk = k(dnib == 1); % dnib = 1 : one 30-bit difference
if ~isempty(kk)
ddd(1,kk) = bitsign(frame32(kk),30);
end
kk = k(dnib == 2); % dnib = 2 : two 15-bit differences
if ~isempty(kk)
ddd(1:2,kk) = bitsplit(frame32(kk),30,15);
end
kk = k(dnib == 3); % dnib = 3 : three 10-bit differences
if ~isempty(kk)
ddd(1:3,kk) = bitsplit(frame32(kk),30,10);
end
end
k = find(nibbles == 3); % nibble = 3 : must look in dnib
if ~isempty(k)
dnib = bitshift(frame32(k),-30);
kk = k(dnib == 0); % dnib = 0 : five 6-bit difference
if ~isempty(kk)
ddd(1:5,kk) = bitsplit(frame32(kk),30,6);
end
kk = k(dnib == 1); % dnib = 1 : six 5-bit differences
if ~isempty(kk)
ddd(1:6,kk) = bitsplit(frame32(kk),30,5);
end
kk = k(dnib == 2); % dnib = 2 : seven 4-bit differences (28 bits!)
if ~isempty(kk)
ddd(1:7,kk) = bitsplit(frame32(kk),28,4);
end
end
case 3 % *** STEIM-3 DECODING IS ALPHA AND UNTESTED ***
% STEIM-3: 7 cases following the nibbles
ddd = NaN*ones(9,numel(frame32)); % initiates array with NaN
k = find(nibbles == 0); % nibble = 0 : two 16-bit differences
if ~isempty(k)
ddd(1:2,k) = bitsplit(frame32(k),32,16);
end
k = find(nibbles == 1); % nibble = 1 : four 8-bit differences
if ~isempty(k)
ddd(1:4,k) = bitsplit(frame32(k),32,8);
end
k = find(nibbles == 2); % nibble = 2 : must look even dnib
if ~isempty(k)
dnib2 = bitshift(frame32(k(2:2:end)),-30);
w60 = bitand(frame32(k(2:2:end)),1073741823) …
+ bitshift(bitand(frame32(k(1:2:end)),1073741823),30); % concatenates two 30-bit words
kk = find(dnib2 == 0); % dnib = 0: five 12-bit differences (60 bits)
if ~isempty(kk)
ddd(1:5,k(2*kk)) = bitsplit(w60,60,12);
end
kk = find(dnib2 == 1); % dnib = 1: three 20-bit differences (60 bits)
if ~isempty(kk)
ddd(1:3,k(2*kk)) = bitsplit(w60,60,20);
end
end
k = find(nibbles == 3); % nibble = 3 : must look 3rd bit
if ~isempty(k)
dnib = bitshift(frame32(k),-27);
kk = k(dnib == 24); % dnib = 11000 : nine 3-bit differences (27 bits)
if ~isempty(kk)
ddd(1:9,kk) = bitsplit(frame32(kk),27,3);
end
kk = k(dnib == 25); % dnib = 11001 : Not A Difference
if ~isempty(kk)
ddd(1,kk) = bitsign(frame32(kk),27);
end
kk = k(dnib > 27); % dnib = 111.. : 29-bit sample (29 bits)
if ~isempty(kk)
ddd(1,kk) = bitsign(frame32(kk),29);
end
end
end
% Little-endian coding: needs to swap bytes
if ~WordOrder
ddd = flipud(ddd);
end
dd = ddd(~isnan(ddd)); % reduces initial array ddd: dd is non-NaN values of ddd
% controls the number of samples
if numel(dd) ~= D.NumberSamples
if verbose > 1
warning(‘RDMSEED:DataIntegrity’,’Problem in %s sequence # %s [%d-%03d %02d:%02d:%07.4f]: number of samples in header (%d) does not equal data (%d).n’, …
D.EncodingFormatName{:},D.SequenceNumber,D.RecordStartTimeISO,D.NumberSamples,numel(dd));
end
if numel(dd) < D.NumberSamples
D.NumberSamples = numel(dd);
end
end
% rebuilds the data vector by integrating the differences
D.d = cumsum([x0;dd(2:D.NumberSamples)]);
% controls data integrity…
if D.d(end) ~= xn
warning(‘RDMSEED:DataIntegrity’,’Problem in %s sequence # %s [%s]: data integrity check failed, last_data=%d, Xn=%d.n’, …
D.EncodingFormatName{:},D.SequenceNumber,D.RecordStartTimeISO,D.d(end),xn);
end
if D.NumberSamples == 0
D.d = nan(0,1);
end
% for debug purpose…
if verbose > 2
D.dd = dd;
D.nibbles = nibbles;
D.x0 = x0;
D.xn = xn;
end
case 12
% — decoding format: GEOSCOPE multiplexed 24-bit integer
D.EncodingFormatName = {‘GEOSCOPE24’};
dd = fread(fid,(D.DataRecordSize – D.OffsetBeginData)/3,’bit24=>double’);
if xor(~WordOrder,le)
dd = swapbytes(dd);
end
D.d = dd(1:D.NumberSamples);
case {13,14}
% — decoding format: GEOSCOPE multiplexed 16/3 and 16/4 bit gain ranged
% (13): 16/3-bit (bit 15 is unused)
% (14): 16/4-bit
% bits 15-12 = 3 or 4-bit gain exponent (positive)
% bits 11-0 = 12-bit mantissa (positive)
% => data = (mantissa – 2048) / 2^gain
geoscope = 7 + 8*(EncodingFormat==14); % mask for gain exponent
D.EncodingFormatName = {sprintf(‘GEOSCOPE16-%d’,EncodingFormat-10)};
dd = fread(fid,(D.DataRecordSize – D.OffsetBeginData)/2,’*uint16′);
if xor(~WordOrder,le)
dd = swapbytes(dd);
end
dd = (double(bitand(dd,2^12-1))-2^11)./2.^double(bitand(bitshift(dd,-12),geoscope));
D.d = dd(1:D.NumberSamples);
case 15
% — decoding format: US National Network compression
D.EncodingFormatName = {‘USNN’};
uncoded = 1;
case 16
% — decoding format: CDSN 16-bit gain ranged
D.EncodingFormatName = {‘CDSN’};
uncoded = 1;
case 17
% — decoding format: Graefenberg 16-bit gain ranged
D.EncodingFormatName = {‘GRAEFENBERG’};
uncoded = 1;
case 18
% — decoding format: IPG – Strasbourg 16-bit gain ranged
D.EncodingFormatName = {‘IPGS’};
uncoded = 1;
case 30
% — decoding format: SRO format
D.EncodingFormatName = {‘SRO’};
uncoded = 1;
case 31
% — decoding format: HGLP format
D.EncodingFormatName = {‘HGLP’};
uncoded = 1;
case 32
% — decoding format: DWWSSN gain ranged format
D.EncodingFormatName = {‘DWWSSN’};
uncoded = 1;
case 33
% — decoding format: RSTN 16-bit gain ranged
D.EncodingFormatName = {‘RSTN’};
uncoded = 1;
otherwise
D.EncodingFormatName = {sprintf(‘** Unknown (%d) **’,EncodingFormat)};
uncoded = 1;
end
if uncoded
error(‘Sorry, the encoding format "%s" is not yet implemented.’,D.EncodingFormatName);
end
% Applies time correction (if needed)
D.RecordStartTimeMATLAB = datenum(double([D.RecordStartTime(1),0,D.RecordStartTime(2:5)])) …
+ (~notc & bitand(D.ActivityFlags,2) == 0)*D.TimeCorrection/1e4/86400;
tv = datevec(D.RecordStartTimeMATLAB);
doy = datenum(tv(1:3)) – datenum(tv(1),1,0);
D.RecordStartTime = [tv(1),doy,tv(4:5),round(tv(6)*1e4)/1e4];
D.RecordStartTimeISO = sprintf(‘%4d-%03d %02d:%02d:%07.4f’,D.RecordStartTime);
D.t = D.RecordStartTimeMATLAB;
% makes the time vector and applies time correction (if needed)
if EncodingFormat > 0
D.t = D.t + (0:(D.NumberSamples-1))’/(D.SampleRate*86400);
end
offset = ftell(fid);
fread(fid,1,’char’); % this is to force EOF=1 on last record.
if feof(fid)
offset = -1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function c = splitfield(s,d)
% splitfield(S) splits string S of D-character separated field names
C = textscan(s,’%s’,’Delimiter’,d);
c = C{1};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [d,swapflag] = readbtime
% readbtime reads BTIME structure from current opened file and returns
% D = [YEAR,DAY,HOUR,MINUTE,SECONDS]
global fid forcebe
Year = fread(fid,1,’*uint16′);
DayOfYear = fread(fid,1,’*uint16′);
Hours = fread(fid,1,’uint8′);
Minutes = fread(fid,1,’uint8′);
Seconds = fread(fid,1,’uint8′);
fseek(fid,1,0); % skip 1 byte (unused)
Seconds0001 = fread(fid,1,’*uint16′);
% Automatic detection of little/big-endian encoding
% — by F. Beauducel, March 2014 —
%
% If the 2-byte day is >= 512, the file is not opened in the correct
% endianness. If the day is 1 or 256, there is a possible byte-swap and we
% need to check also the year; but we need to consider what is a valid year:
% – years from 1801 to 2047 are OK (swapbytes >= 2312)
% – years from 2048 to 2055 are OK (swapbytes <= 1800)
% – year 2056 is ambiguous (swapbytes = 2056)
% – years from 2057 to 2311 are OK (swapbytes >= 2312)
% – year 1799 is ambiguous (swapbytes = 1799)
% – year 1800 is suspicious (swapbytes = 2055)
%
% Thus, the only cases for which we are ‘sure’ there is a byte-swap, are:
% – day >= 512
% – (day == 1 or day == 256) and (year < 1799 or year > 2311)
%
% Note: in IRIS libmseed, the test is only year>2050 or year<1920.
if ~forcebe && (DayOfYear >= 512 || (ismember(DayOfYear,[1,256]) && (Year > 2311 || Year < 1799)))
swapflag = 1;
Year = swapbytes(Year);
DayOfYear = swapbytes(DayOfYear);
Seconds0001 = swapbytes(Seconds0001);
else
swapflag = 0;
end
d = [double(Year),double(DayOfYear),Hours,Minutes,Seconds + double(Seconds0001)/1e4];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function d = bitsplit(x,b,n)
% bitsplit(X,B,N) splits the B-bit number X into signed N-bit array
% X must be unsigned integer class
% N ranges from 1 to B
% B is a multiple of N
sign = repmat((b:-n:n)’,1,size(x,1));
x = repmat(x’,b/n,1);
d = double(bitand(bitshift(x,flipud(sign-b)),2^n-1)) …
– double(bitget(x,sign))*2^n;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function d = bitsign(x,n)
% bitsign(X,N) returns signed double value from unsigned N-bit number X.
% This is equivalent to bitsplit(X,N,N), but the formula is simplified so
% it is much more efficient
d = double(bitand(x,2^n-1)) – double(bitget(x,n)).*2^n; How can i change the subplots sequence in rdmseed.m, while plotting a miniseed file with e.g. 5 stations?
function varargout = rdmseed(varargin)
if nargin > 6
error(‘Too many input arguments.’)
end
% global variables shared with sub-functions
global f fid offset le ef wo rl forcebe verbose notc force
% default input arguments
makeplot = 0; % make plot flag
verbose = 0; % verbose flag/level
forcebe = 0; % force big-endian
ef = 10; % encoding format default
wo = 1; % word order default
rl = 2^12; % record length default
force = 0; % force input argument over blockette 1000 (UNDOCUMENTED)
notc = 0; % force no time correction (over ActivityFlags)
nullhead = 0; % allow null bytes before header
if nargin < 1
[filename,pathname] = uigetfile(‘*’,’Please select a miniSEED file…’);
f = fullfile(pathname,filename);
else
f = varargin{1};
end
if ~ischar(f) || ~exist(f,’file’)
error(‘File %s does not exist.’,f);
end
if nargin > 1
verbose = any(strcmpi(varargin,’v’)) + 2*any(strcmpi(varargin,’vv’)) …
+ 3*any(strcmpi(varargin,’vvv’));
makeplot = any(strcmpi(varargin,’plot’));
forcebe = any(strcmpi(varargin,’be’));
notc = any(strcmpi(varargin,’notc’));
force = any(strcmpi(varargin,’force’));
nullhead = any(strcmpi(varargin,’nullhead’));
end
nargs = (makeplot>0) + (verbose>0) + (forcebe>0) + (notc>0) + (force>0) …
+ (nullhead>0);
if nargin > (1 + nargs)
ef = varargin{2};
if ~isnumeric(ef) || ~any(ef==[0:5,10:19,30:33])
error(‘Argument ENCODINGFORMAT must be a valid FDSN code value.’);
end
end
if nargin > (2 + nargs)
wo = varargin{3};
if ~isnumeric(wo) || (wo ~= 0 && wo ~= 1)
error(‘Argument WORDORDER must be 0 or 1.’);
end
end
if nargin > (3 + nargs)
rl = varargin{4};
if ~isnumeric(rl) || rl < 256 || rem(log(rl)/log(2),1) ~= 0
error(‘Argument RECORDLENGTH must be a power of 2 and greater or equal to 256.’);
end
end
if nargout == 0
makeplot = 1;
end
% sensible limits for multiplexed files
max_channels = 20; % absolute max number of channels to plot
max_channel_label = 6; % max. number of channels for y-labels
% file is opened in Big-Endian encoding (this is encouraged by SEED)
fid = fopen(f,’rb’,’ieee-be’);
le = 0;
offset = 0;
% — tests if the header is mini-SEED
% the 7th character must be one of the "data header/quality indicator", usually ‘D’
header = fread(fid,20,’*char’);
if ~ismember(header(7),’DRMQ’)
if ismember(header(7),’VAST’)
error(‘File seems to be a SEED Volume. Cannot read it.’);
else
if header(1)==0
if nullhead
if verbose
fprintf(‘Null header option: bypassing…’);
end
c = 0;
fseek(fid,0,’bof’);
while c==0
c = fread(fid,1,’*char’);
offset = offset + 1;
end
if verbose
fprintf(‘ %d null bytes.n’,offset);
end
header = fread(fid,6,’*char’);
if ~ismember(header(6),’DRMQ’)
error(‘File is not in mini-SEED format. Cannot read it.’);
else
offset = offset – 1;
end
else
error(‘File starts with null bytes… if you believe it is still a miniseed file, try the ”nullhead” option.’);
end
else
error(‘File is not in mini-SEED format. Cannot read it.’);
end
end
end
i = 1;
while offset >= 0
X(i) = read_data_record;
i = i + 1;
end
fclose(fid);
if nargout > 0
varargout{1} = X;
end
% — analyses data
if makeplot || nargout > 1
% test if the file is multiplexed or a single channel
un = unique(cellstr(char(X.ChannelFullName)));
nc = numel(un);
for i = 1:nc
k = find(strcmp(cellstr(char(X.ChannelFullName)),un{i}));
I(i).ChannelFullName = X(k(1)).ChannelFullName;
I(i).XBlockIndex = k;
I(i).ClockDrift = ([diff(cat(1,X(k).RecordStartTimeMATLAB));NaN]*86400 – cat(1,X(k).NumberSamples)./cat(1,X(k).SampleRate))./cat(1,X(k).NumberSamples);
I(i).OverlapBlockIndex = k(find(I(i).ClockDrift.*cat(1,X(k).NumberSamples).*cat(1,X(k).SampleRate) < -.5) + 1);
I(i).OverlapTime = cat(1,X(I(i).OverlapBlockIndex).RecordStartTimeMATLAB);
I(i).GapBlockIndex = k(find(I(i).ClockDrift.*cat(1,X(k).NumberSamples).*cat(1,X(k).SampleRate) > .5) + 1);
I(i).GapTime = cat(1,X(I(i).GapBlockIndex).RecordStartTimeMATLAB);
end
end
if nargout > 1
varargout{2} = I;
end
% — plots the data
if makeplot
figure
xlim = [min(cat(1,X.t)),max(cat(1,X.t))];
% test if all data records have the same length
rl = unique(cat(1,X.DataRecordSize));
if numel(rl) == 1
rl_text = sprintf(‘%d bytes’,rl);
else
rl_text = sprintf(‘%d-%d bytes’,min(rl),max(rl));
end
% test if all data records have the same sampling rate
sr = unique(cat(1,X.SampleRate));
if numel(sr) == 1
sr_text = sprintf(‘%g Hz’,sr);
else
sr_text = sprintf(‘%d # samp. rates’,numel(sr));
end
% test if all data records have the same encoding format
ef = unique(cellstr(cat(1,X.EncodingFormatName)));
if numel(ef) == 1
ef_text = sprintf(‘%s’,ef{:});
else
ef_text = sprintf(‘%d different encod. formats’,numel(ef));
end
if nc == 1
plot(cat(1,X.t),cat(1,X.d))
hold on
for i = 1:length(I.GapBlockIndex)
plot(I.GapTime(i),X(I.GapBlockIndex(i)).d(1),’*r’)
end
for i = 1:length(I.OverlapBlockIndex)
plot(I.OverlapTime(i),X(I.OverlapBlockIndex(i)).d(1),’og’)
end
hold off
set(gca,’XLim’,xlim)
datetick(‘x’,’keeplimits’)
grid on
xlabel(sprintf(‘Timen(%s to %s)’,datestr(xlim(1)),datestr(xlim(2))))
ylabel(‘Counts’)
title(sprintf(‘mini-SEED file "%s"n%s (%d rec. @ %s – %g samp. @ %s – %s)’, …
f,un{1},length(X),rl_text,numel(cat(1,X.d)),sr_text,ef_text),’Interpreter’,’none’)
else
% plot is done only for real data channels…
if nc > max_channels
warning(‘Plot has been limited to %d channels (over %d). See help to manage multiplexed file.’, …
max_channels,nc);
nc = max_channels;
end
for i = 1:nc
subplot(nc*2,1,i*2 + (-1:0))
k = I(i).XBlockIndex;
if ~any(strcmp(‘ASCII’,cellstr(cat(1,X(k).EncodingFormatName))))
plot(cat(1,X(k).t),cat(1,X(k).d))
hold on
for ii = 1:length(I(i).GapBlockIndex)
if ~isempty(X(I(i).GapBlockIndex(ii)).d)
plot(I(i).GapTime(ii),X(I(i).GapBlockIndex(ii)).d,’r’)
else
plot(repmat(I(i).GapTime(ii),1,2),ylim,’r’)
end
end
for ii = 1:length(I(i).OverlapBlockIndex)
if ~isempty(X(I(i).OverlapBlockIndex(ii)).d)
plot(I(i).OverlapTime(ii),X(I(i).OverlapBlockIndex(ii)).d,’g’)
else
plot(repmat(I(i).OverlapTime(ii),1,2),ylim,’g’)
end
end
hold off
end
set(gca,’XLim’,xlim,’FontSize’,8)
h = ylabel(un{i},’Interpreter’,’none’);
if nc > max_channel_label
set(gca,’YTick’,[])
set(h,’Rotation’,0,’HorizontalAlignment’,’right’,’FontSize’,8)
end
datetick(‘x’,’keeplimits’)
set(gca,’XTickLabel’,[])
grid on
if i == 1
title(sprintf(‘mini-SEED file "%s"n%d channels (%d rec. @ %s – %g data – %s – %s)’, …
f,length(un),length(X),rl_text,numel(cat(1,X(k).d)),sr_text,ef_text),’Interpreter’,’none’)
end
if i == nc
datetick(‘x’,’keeplimits’)
xlabel(sprintf(‘Timen(%s to %s)’,datestr(xlim(1)),datestr(xlim(2))))
end
end
v = version;
if str2double(v(1))>=7
linkaxes(findobj(gcf,’type’,’axes’),’x’)
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function D = read_data_record
% read_data_record uses global variables f, fid, offset, le, ef, wo, rl,
% and verbose. It reads a data record and returns a structure D.
global f fid offset le ef wo rl verbose notc force
fseek(fid,offset,’bof’);
% — read fixed section of Data Header (48 bytes)
D.SequenceNumber = fread(fid,6,’*char’)’;
D.DataQualityIndicator = fread(fid,1,’*char’);
D.ReservedByte = fread(fid,1,’*char’);
D.StationIdentifierCode = fread(fid,5,’*char’)’;
D.LocationIdentifier = fread(fid,2,’*char’)’;
D.ChannelIdentifier = fread(fid,3,’*char’)’;
D.NetworkCode = fread(fid,2,’*char’)’;
D.ChannelFullName = sprintf(‘%s:%s:%s:%s’,deblank(D.NetworkCode), …
deblank(D.StationIdentifierCode),deblank(D.LocationIdentifier), …
deblank(D.ChannelIdentifier));
% Start Time decoding
[D.RecordStartTime,swapflag] = readbtime;
D.RecordStartTimeISO = sprintf(‘%4d-%03d %02d:%02d:%07.4f’,D.RecordStartTime);
if swapflag
if le
machinefmt = ‘ieee-be’;
le = 0;
else
machinefmt = ‘ieee-le’;
le = 1;
end
position = ftell(fid);
fclose(fid);
fid = fopen(f,’rb’,machinefmt);
fseek(fid,position,’bof’);
if verbose > 0
warning(‘RDMSEED:DataIntegrity’, …
‘Sequence # %s: need to switch file encoding to %s…n’, …
D.SequenceNumber,machinefmt);
end
end
D.NumberSamples = fread(fid,1,’uint16′);
% Sample Rate decoding
SampleRateFactor = fread(fid,1,’int16′);
SampleRateMultiplier = fread(fid,1,’int16′);
if SampleRateFactor > 0
if SampleRateMultiplier >= 0
D.SampleRate = SampleRateFactor*SampleRateMultiplier;
else
D.SampleRate = -1*SampleRateFactor/SampleRateMultiplier;
end
else
if SampleRateMultiplier >= 0
D.SampleRate = -1*SampleRateMultiplier/SampleRateFactor;
else
D.SampleRate = 1/(SampleRateFactor*SampleRateMultiplier);
end
end
D.ActivityFlags = fread(fid,1,’uint8′);
D.IOFlags = fread(fid,1,’uint8′);
D.DataQualityFlags = fread(fid,1,’uint8′);
D.NumberBlockettesFollow = fread(fid,1,’uint8′);
D.TimeCorrection = fread(fid,1,’int32′); % Time correction in 0.0001 s
D.OffsetBeginData = fread(fid,1,’uint16′);
D.OffsetFirstBlockette = fread(fid,1,’uint16′);
% — read the blockettes
OffsetNextBlockette = D.OffsetFirstBlockette;
D.BLOCKETTES = [];
b2000 = 0; % Number of Blockette 2000
for i = 1:D.NumberBlockettesFollow
fseek(fid,offset + OffsetNextBlockette,’bof’);
BlocketteType = fread(fid,1,’uint16′);
switch BlocketteType
case 1000
% BLOCKETTE 1000 = Data Only SEED (8 bytes)
OffsetNextBlockette = fread(fid,1,’uint16′);
D.BLOCKETTES.B1000.EncodingFormat = fread(fid,1,’uint8′);
D.BLOCKETTES.B1000.WordOrder = fread(fid,1,’uint8′);
D.BLOCKETTES.B1000.DataRecordLength = fread(fid,1,’uint8′);
D.BLOCKETTES.B1000.Reserved = fread(fid,1,’uint8′);
case 1001
% BLOCKETTE 1001 = Data Extension (8 bytes)
OffsetNextBlockette = fread(fid,1,’uint16′);
D.BLOCKETTES.B1001.TimingQuality = fread(fid,1,’uint8′);
D.BLOCKETTES.B1001.Micro_sec = fread(fid,1,’int8′);
D.BLOCKETTES.B1001.Reserved = fread(fid,1,’uint8′);
D.BLOCKETTES.B1001.FrameCount = fread(fid,1,’uint8′);
case 100
% BLOCKETTE 100 = Sample Rate (12 bytes)
OffsetNextBlockette = fread(fid,1,’uint16′);
D.BLOCKETTES.B100.ActualSampleRate = fread(fid,1,’float32′);
D.BLOCKETTES.B100.Flags = fread(fid,1,’uint8′);
D.BLOCKETTES.B100.Reserved = fread(fid,1,’uint8′);
case 500
% BLOCKETTE 500 = Timing (200 bytes)
OffsetNextBlockette = fread(fid,1,’uint16′);
D.BLOCKETTES.B500.VCOCorrection = fread(fid,1,’float32′);
D.BLOCKETTES.B500.TimeOfException = readbtime;
D.BLOCKETTES.B500.MicroSec = fread(fid,1,’int8′);
D.BLOCKETTES.B500.ReceptionQuality = fread(fid,1,’uint8′);
D.BLOCKETTES.B500.ExceptionCount = fread(fid,1,’uint16′);
D.BLOCKETTES.B500.ExceptionType = fread(fid,16,’*char’)’;
D.BLOCKETTES.B500.ClockModel = fread(fid,32,’*char’)’;
D.BLOCKETTES.B500.ClockStatus = fread(fid,128,’*char’)’;
case 2000
% BLOCKETTE 2000 = Opaque Data (variable length)
b2000 = b2000 + 1;
OffsetNextBlockette = fread(fid,1,’uint16′);
BlocketteLength = fread(fid,1,’uint16′);
OffsetOpaqueData = fread(fid,1,’uint16′);
D.BLOCKETTES.B2000(b2000).RecordNumber = fread(fid,1,’uint32′);
D.BLOCKETTES.B2000(b2000).DataWordOrder = fread(fid,1,’uint8′);
D.BLOCKETTES.B2000(b2000).Flags = fread(fid,1,’uint8′);
NumberHeaderFields = fread(fid,1,’uint8′);
HeaderFields = splitfield(fread(fid,OffsetOpaqueData-15,’*char’)’,’~’);
D.BLOCKETTES.B2000(b2000).HeaderFields = HeaderFields(1:NumberHeaderFields);
% Opaque data are stored as a single char string, but must be
% decoded using appropriate format (e.g., Quanterra Q330)
D.BLOCKETTES.B2000(b2000).OpaqueData = fread(fid,BlocketteLength-OffsetOpaqueData,’*char’)’;
otherwise
OffsetNextBlockette = fread(fid,1,’uint16′);
if verbose > 0
warning(‘RDMSEED:UnknownBlockette’, …
‘Unknown Blockette number %d (%s)!n’, …
BlocketteType,D.ChannelFullName);
end
end
end
% — read the data stream
fseek(fid,offset + D.OffsetBeginData,’bof’);
if ~force && isfield(D.BLOCKETTES,’B1000′)
EncodingFormat = D.BLOCKETTES.B1000.EncodingFormat;
WordOrder = D.BLOCKETTES.B1000.WordOrder;
D.DataRecordSize = 2^D.BLOCKETTES.B1000.DataRecordLength;
else
EncodingFormat = ef;
WordOrder = wo;
D.DataRecordSize = rl;
end
uncoded = 0;
D.d = NaN;
D.t = NaN;
switch EncodingFormat
case 0
% — decoding format: ASCII text
D.EncodingFormatName = {‘ASCII’};
D.d = fread(fid,D.DataRecordSize – D.OffsetBeginData,’*char’)’;
case 1
% — decoding format: 16-bit integers
D.EncodingFormatName = {‘INT16’};
dd = fread(fid,ceil((D.DataRecordSize – D.OffsetBeginData)/2),’*int16′);
if xor(~WordOrder,le)
dd = swapbytes(dd);
end
D.d = dd(1:D.NumberSamples);
case 2
% — decoding format: 24-bit integers
D.EncodingFormatName = {‘INT24’};
dd = fread(fid,ceil((D.DataRecordSize – D.OffsetBeginData)/3),’bit24=>int32′);
if xor(~WordOrder,le)
dd = swapbytes(dd);
end
D.d = dd(1:D.NumberSamples);
case 3
% — decoding format: 32-bit integers
D.EncodingFormatName = {‘INT32’};
dd = fread(fid,ceil((D.DataRecordSize – D.OffsetBeginData)/4),’*int32′);
if xor(~WordOrder,le)
dd = swapbytes(dd);
end
D.d = dd(1:D.NumberSamples);
case 4
% — decoding format: IEEE floating point
D.EncodingFormatName = {‘FLOAT32’};
dd = fread(fid,ceil((D.DataRecordSize – D.OffsetBeginData)/4),’*float’);
if xor(~WordOrder,le)
dd = swapbytes(dd);
end
D.d = dd(1:D.NumberSamples);
case 5
% — decoding format: IEEE double precision floating point
D.EncodingFormatName = {‘FLOAT64’};
dd = fread(fid,ceil((D.DataRecordSize – D.OffsetBeginData)/8),’*double’);
if xor(~WordOrder,le)
dd = swapbytes(dd);
end
D.d = dd(1:D.NumberSamples);
case {10,11,19}
% — decoding formats: STEIM-1 and STEIM-2 compression
% (c) Joseph M. Steim, Quanterra Inc., 1994
steim = find(EncodingFormat==[10,11,19]);
D.EncodingFormatName = {sprintf(‘STEIM%d’,steim)};
% Steim compression decoding strategy optimized for Matlab
% — by F. Beauducel, October 2010 —
%
% 1. loads all data into a single 16xM uint32 array
% 2. gets all nibbles from the first row splitted into 2-bit values
% 3. for each possible nibble value, selects (find) and decodes
% (bitsplit) all the corresponding words, and stores results
% in a 4xN (STEIM1) or 7xN (STEIM2) array previously filled with
% NaN’s. For STEIM2 with nibbles 2 or 3, decodes also dnib values
% (first 2-bit of the word)
% 5. reduces this array with non-NaN values only
% 6. integrates with cumsum
%
% This method is about 30 times faster than a ‘C-like’ loops coding…
frame32 = fread(fid,[16,(D.DataRecordSize – D.OffsetBeginData)/64],’*uint32′);
if xor(~WordOrder,le)
frame32 = swapbytes(frame32);
end
% specific processes for STEIM-3
if steim == 3
% first bit = 1 means second differences
SecondDiff = bitshift(frame32(1,:),-31);
% checks for "squeezed flag"… and replaces frame32(1,:)
squeezed = bitand(bitshift(frame32(1,:),-24),127);
k = find(bitget(squeezed,7));
if ~isempty(k)
moredata24 = bitand(frame32(1,k),16777215);
k = find(squeezed == 80); % upper nibble 8-bit = 0x50
if ~isempty(k)
frame32(1,k) = hex2dec(‘15555555’);
end
k = find(squeezed == 96); % upper nibble 8-bit = 0x60
if ~isempty(k)
frame32(1,k) = hex2dec(‘2aaaaaaa’);
end
k = find(squeezed == 112); % upper nibble 8-bit = 0x70
if ~isempty(k)
frame32(1,k) = hex2dec(‘3fffffff’);
end
end
end
% nibbles is an array of the same size as frame32…
nibbles = bitand(bitshift(repmat(frame32(1,:),16,1),repmat(-30:2:0,size(frame32,2),1)’),3);
x0 = bitsign(frame32(2,1),32); % forward integration constant
xn = bitsign(frame32(3,1),32); % reverse integration constant
switch steim
case 1
% STEIM-1: 3 cases following the nibbles
ddd = NaN*ones(4,numel(frame32)); % initiates array with NaN
k = find(nibbles == 1); % nibble = 1 : four 8-bit differences
if ~isempty(k)
ddd(1:4,k) = bitsplit(frame32(k),32,8);
end
k = find(nibbles == 2); % nibble = 2 : two 16-bit differences
if ~isempty(k)
ddd(1:2,k) = bitsplit(frame32(k),32,16);
end
k = find(nibbles == 3); % nibble = 3 : one 32-bit difference
if ~isempty(k)
ddd(1,k) = bitsign(frame32(k),32);
end
case 2
% STEIM-2: 7 cases following the nibbles and dnib
ddd = NaN*ones(7,numel(frame32)); % initiates array with NaN
k = find(nibbles == 1); % nibble = 1 : four 8-bit differences
if ~isempty(k)
ddd(1:4,k) = bitsplit(frame32(k),32,8);
end
k = find(nibbles == 2); % nibble = 2 : must look in dnib
if ~isempty(k)
dnib = bitshift(frame32(k),-30);
kk = k(dnib == 1); % dnib = 1 : one 30-bit difference
if ~isempty(kk)
ddd(1,kk) = bitsign(frame32(kk),30);
end
kk = k(dnib == 2); % dnib = 2 : two 15-bit differences
if ~isempty(kk)
ddd(1:2,kk) = bitsplit(frame32(kk),30,15);
end
kk = k(dnib == 3); % dnib = 3 : three 10-bit differences
if ~isempty(kk)
ddd(1:3,kk) = bitsplit(frame32(kk),30,10);
end
end
k = find(nibbles == 3); % nibble = 3 : must look in dnib
if ~isempty(k)
dnib = bitshift(frame32(k),-30);
kk = k(dnib == 0); % dnib = 0 : five 6-bit difference
if ~isempty(kk)
ddd(1:5,kk) = bitsplit(frame32(kk),30,6);
end
kk = k(dnib == 1); % dnib = 1 : six 5-bit differences
if ~isempty(kk)
ddd(1:6,kk) = bitsplit(frame32(kk),30,5);
end
kk = k(dnib == 2); % dnib = 2 : seven 4-bit differences (28 bits!)
if ~isempty(kk)
ddd(1:7,kk) = bitsplit(frame32(kk),28,4);
end
end
case 3 % *** STEIM-3 DECODING IS ALPHA AND UNTESTED ***
% STEIM-3: 7 cases following the nibbles
ddd = NaN*ones(9,numel(frame32)); % initiates array with NaN
k = find(nibbles == 0); % nibble = 0 : two 16-bit differences
if ~isempty(k)
ddd(1:2,k) = bitsplit(frame32(k),32,16);
end
k = find(nibbles == 1); % nibble = 1 : four 8-bit differences
if ~isempty(k)
ddd(1:4,k) = bitsplit(frame32(k),32,8);
end
k = find(nibbles == 2); % nibble = 2 : must look even dnib
if ~isempty(k)
dnib2 = bitshift(frame32(k(2:2:end)),-30);
w60 = bitand(frame32(k(2:2:end)),1073741823) …
+ bitshift(bitand(frame32(k(1:2:end)),1073741823),30); % concatenates two 30-bit words
kk = find(dnib2 == 0); % dnib = 0: five 12-bit differences (60 bits)
if ~isempty(kk)
ddd(1:5,k(2*kk)) = bitsplit(w60,60,12);
end
kk = find(dnib2 == 1); % dnib = 1: three 20-bit differences (60 bits)
if ~isempty(kk)
ddd(1:3,k(2*kk)) = bitsplit(w60,60,20);
end
end
k = find(nibbles == 3); % nibble = 3 : must look 3rd bit
if ~isempty(k)
dnib = bitshift(frame32(k),-27);
kk = k(dnib == 24); % dnib = 11000 : nine 3-bit differences (27 bits)
if ~isempty(kk)
ddd(1:9,kk) = bitsplit(frame32(kk),27,3);
end
kk = k(dnib == 25); % dnib = 11001 : Not A Difference
if ~isempty(kk)
ddd(1,kk) = bitsign(frame32(kk),27);
end
kk = k(dnib > 27); % dnib = 111.. : 29-bit sample (29 bits)
if ~isempty(kk)
ddd(1,kk) = bitsign(frame32(kk),29);
end
end
end
% Little-endian coding: needs to swap bytes
if ~WordOrder
ddd = flipud(ddd);
end
dd = ddd(~isnan(ddd)); % reduces initial array ddd: dd is non-NaN values of ddd
% controls the number of samples
if numel(dd) ~= D.NumberSamples
if verbose > 1
warning(‘RDMSEED:DataIntegrity’,’Problem in %s sequence # %s [%d-%03d %02d:%02d:%07.4f]: number of samples in header (%d) does not equal data (%d).n’, …
D.EncodingFormatName{:},D.SequenceNumber,D.RecordStartTimeISO,D.NumberSamples,numel(dd));
end
if numel(dd) < D.NumberSamples
D.NumberSamples = numel(dd);
end
end
% rebuilds the data vector by integrating the differences
D.d = cumsum([x0;dd(2:D.NumberSamples)]);
% controls data integrity…
if D.d(end) ~= xn
warning(‘RDMSEED:DataIntegrity’,’Problem in %s sequence # %s [%s]: data integrity check failed, last_data=%d, Xn=%d.n’, …
D.EncodingFormatName{:},D.SequenceNumber,D.RecordStartTimeISO,D.d(end),xn);
end
if D.NumberSamples == 0
D.d = nan(0,1);
end
% for debug purpose…
if verbose > 2
D.dd = dd;
D.nibbles = nibbles;
D.x0 = x0;
D.xn = xn;
end
case 12
% — decoding format: GEOSCOPE multiplexed 24-bit integer
D.EncodingFormatName = {‘GEOSCOPE24’};
dd = fread(fid,(D.DataRecordSize – D.OffsetBeginData)/3,’bit24=>double’);
if xor(~WordOrder,le)
dd = swapbytes(dd);
end
D.d = dd(1:D.NumberSamples);
case {13,14}
% — decoding format: GEOSCOPE multiplexed 16/3 and 16/4 bit gain ranged
% (13): 16/3-bit (bit 15 is unused)
% (14): 16/4-bit
% bits 15-12 = 3 or 4-bit gain exponent (positive)
% bits 11-0 = 12-bit mantissa (positive)
% => data = (mantissa – 2048) / 2^gain
geoscope = 7 + 8*(EncodingFormat==14); % mask for gain exponent
D.EncodingFormatName = {sprintf(‘GEOSCOPE16-%d’,EncodingFormat-10)};
dd = fread(fid,(D.DataRecordSize – D.OffsetBeginData)/2,’*uint16′);
if xor(~WordOrder,le)
dd = swapbytes(dd);
end
dd = (double(bitand(dd,2^12-1))-2^11)./2.^double(bitand(bitshift(dd,-12),geoscope));
D.d = dd(1:D.NumberSamples);
case 15
% — decoding format: US National Network compression
D.EncodingFormatName = {‘USNN’};
uncoded = 1;
case 16
% — decoding format: CDSN 16-bit gain ranged
D.EncodingFormatName = {‘CDSN’};
uncoded = 1;
case 17
% — decoding format: Graefenberg 16-bit gain ranged
D.EncodingFormatName = {‘GRAEFENBERG’};
uncoded = 1;
case 18
% — decoding format: IPG – Strasbourg 16-bit gain ranged
D.EncodingFormatName = {‘IPGS’};
uncoded = 1;
case 30
% — decoding format: SRO format
D.EncodingFormatName = {‘SRO’};
uncoded = 1;
case 31
% — decoding format: HGLP format
D.EncodingFormatName = {‘HGLP’};
uncoded = 1;
case 32
% — decoding format: DWWSSN gain ranged format
D.EncodingFormatName = {‘DWWSSN’};
uncoded = 1;
case 33
% — decoding format: RSTN 16-bit gain ranged
D.EncodingFormatName = {‘RSTN’};
uncoded = 1;
otherwise
D.EncodingFormatName = {sprintf(‘** Unknown (%d) **’,EncodingFormat)};
uncoded = 1;
end
if uncoded
error(‘Sorry, the encoding format "%s" is not yet implemented.’,D.EncodingFormatName);
end
% Applies time correction (if needed)
D.RecordStartTimeMATLAB = datenum(double([D.RecordStartTime(1),0,D.RecordStartTime(2:5)])) …
+ (~notc & bitand(D.ActivityFlags,2) == 0)*D.TimeCorrection/1e4/86400;
tv = datevec(D.RecordStartTimeMATLAB);
doy = datenum(tv(1:3)) – datenum(tv(1),1,0);
D.RecordStartTime = [tv(1),doy,tv(4:5),round(tv(6)*1e4)/1e4];
D.RecordStartTimeISO = sprintf(‘%4d-%03d %02d:%02d:%07.4f’,D.RecordStartTime);
D.t = D.RecordStartTimeMATLAB;
% makes the time vector and applies time correction (if needed)
if EncodingFormat > 0
D.t = D.t + (0:(D.NumberSamples-1))’/(D.SampleRate*86400);
end
offset = ftell(fid);
fread(fid,1,’char’); % this is to force EOF=1 on last record.
if feof(fid)
offset = -1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function c = splitfield(s,d)
% splitfield(S) splits string S of D-character separated field names
C = textscan(s,’%s’,’Delimiter’,d);
c = C{1};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [d,swapflag] = readbtime
% readbtime reads BTIME structure from current opened file and returns
% D = [YEAR,DAY,HOUR,MINUTE,SECONDS]
global fid forcebe
Year = fread(fid,1,’*uint16′);
DayOfYear = fread(fid,1,’*uint16′);
Hours = fread(fid,1,’uint8′);
Minutes = fread(fid,1,’uint8′);
Seconds = fread(fid,1,’uint8′);
fseek(fid,1,0); % skip 1 byte (unused)
Seconds0001 = fread(fid,1,’*uint16′);
% Automatic detection of little/big-endian encoding
% — by F. Beauducel, March 2014 —
%
% If the 2-byte day is >= 512, the file is not opened in the correct
% endianness. If the day is 1 or 256, there is a possible byte-swap and we
% need to check also the year; but we need to consider what is a valid year:
% – years from 1801 to 2047 are OK (swapbytes >= 2312)
% – years from 2048 to 2055 are OK (swapbytes <= 1800)
% – year 2056 is ambiguous (swapbytes = 2056)
% – years from 2057 to 2311 are OK (swapbytes >= 2312)
% – year 1799 is ambiguous (swapbytes = 1799)
% – year 1800 is suspicious (swapbytes = 2055)
%
% Thus, the only cases for which we are ‘sure’ there is a byte-swap, are:
% – day >= 512
% – (day == 1 or day == 256) and (year < 1799 or year > 2311)
%
% Note: in IRIS libmseed, the test is only year>2050 or year<1920.
if ~forcebe && (DayOfYear >= 512 || (ismember(DayOfYear,[1,256]) && (Year > 2311 || Year < 1799)))
swapflag = 1;
Year = swapbytes(Year);
DayOfYear = swapbytes(DayOfYear);
Seconds0001 = swapbytes(Seconds0001);
else
swapflag = 0;
end
d = [double(Year),double(DayOfYear),Hours,Minutes,Seconds + double(Seconds0001)/1e4];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function d = bitsplit(x,b,n)
% bitsplit(X,B,N) splits the B-bit number X into signed N-bit array
% X must be unsigned integer class
% N ranges from 1 to B
% B is a multiple of N
sign = repmat((b:-n:n)’,1,size(x,1));
x = repmat(x’,b/n,1);
d = double(bitand(bitshift(x,flipud(sign-b)),2^n-1)) …
– double(bitget(x,sign))*2^n;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function d = bitsign(x,n)
% bitsign(X,N) returns signed double value from unsigned N-bit number X.
% This is equivalent to bitsplit(X,N,N), but the formula is simplified so
% it is much more efficient
d = double(bitand(x,2^n-1)) – double(bitget(x,n)).*2^n; rdmseed MATLAB Answers — New Questions
duration of free trial
How long does a free trial of MATLAB last?How long does a free trial of MATLAB last? How long does a free trial of MATLAB last? trial, free trial MATLAB Answers — New Questions
Radar IQ Data FFT/STFT Spectogram Problem
Hello,
I have data which is taken from FMCW radar. Radar specs are as in below:
Center Frequency 77GHz
Bandwidth 160MHz
Range Resolution 1m
PRF 17kHz
Azimuth Beamwidth 1
Scan Rate 10Hz,
mechanical Field of View 9
Output Power 10mW
I use FFT and STFT functions for spectogram plotting. However, it’s seems wrong for me. I want to ask about the usage of those functions. My usage of those functions as in below;
STFT:
window_size = 128;
overlap_length = 120;
fft_points= 256;
window = hamming(window_size);
[S, F, T] = stft(IQ, 1280, ‘Window’, window, ‘OverlapLength’, overlap_length, ‘FFTLength’, fft_points);
new_data_Set_STFT{index_onearray,2} = F; %Frequency
% Amplitude of Spectogram
new_data_Set_STFT{index_onearray,3} =mean(abs(S’));
% IQ data phase
new_data_Set_STFT{index_onearray,4} = atan(imag(S)/real(S));
FFT:
% IQ DATA FFT
new_data_Set{index_onearray,2} = fft(IQ(:,index_global));
% Frequency length
N = length(IQ(:,index_global)); % FFT length
fs = 256; % sampling frequency
calculation = (0:N-1)*(fs/N);
f1 = reshape(calculation, 1280, …
1);
new_data_Set{index_onearray,3} = f1; % Frequency vector
% Amplitude of Spectogram
new_data_Set{index_onearray,4} = abs(IQ(:,index_global));
% IQ data phase
new_data_Set{index_onearray,5} = angle(IQ(:,index_global));
After Plotting, STFT spektogram looks like this. In frequency and time(x time, y frequency) scale it looks wrong. What’s wrong with my code?
FFT:
Thanks for help.Hello,
I have data which is taken from FMCW radar. Radar specs are as in below:
Center Frequency 77GHz
Bandwidth 160MHz
Range Resolution 1m
PRF 17kHz
Azimuth Beamwidth 1
Scan Rate 10Hz,
mechanical Field of View 9
Output Power 10mW
I use FFT and STFT functions for spectogram plotting. However, it’s seems wrong for me. I want to ask about the usage of those functions. My usage of those functions as in below;
STFT:
window_size = 128;
overlap_length = 120;
fft_points= 256;
window = hamming(window_size);
[S, F, T] = stft(IQ, 1280, ‘Window’, window, ‘OverlapLength’, overlap_length, ‘FFTLength’, fft_points);
new_data_Set_STFT{index_onearray,2} = F; %Frequency
% Amplitude of Spectogram
new_data_Set_STFT{index_onearray,3} =mean(abs(S’));
% IQ data phase
new_data_Set_STFT{index_onearray,4} = atan(imag(S)/real(S));
FFT:
% IQ DATA FFT
new_data_Set{index_onearray,2} = fft(IQ(:,index_global));
% Frequency length
N = length(IQ(:,index_global)); % FFT length
fs = 256; % sampling frequency
calculation = (0:N-1)*(fs/N);
f1 = reshape(calculation, 1280, …
1);
new_data_Set{index_onearray,3} = f1; % Frequency vector
% Amplitude of Spectogram
new_data_Set{index_onearray,4} = abs(IQ(:,index_global));
% IQ data phase
new_data_Set{index_onearray,5} = angle(IQ(:,index_global));
After Plotting, STFT spektogram looks like this. In frequency and time(x time, y frequency) scale it looks wrong. What’s wrong with my code?
FFT:
Thanks for help. Hello,
I have data which is taken from FMCW radar. Radar specs are as in below:
Center Frequency 77GHz
Bandwidth 160MHz
Range Resolution 1m
PRF 17kHz
Azimuth Beamwidth 1
Scan Rate 10Hz,
mechanical Field of View 9
Output Power 10mW
I use FFT and STFT functions for spectogram plotting. However, it’s seems wrong for me. I want to ask about the usage of those functions. My usage of those functions as in below;
STFT:
window_size = 128;
overlap_length = 120;
fft_points= 256;
window = hamming(window_size);
[S, F, T] = stft(IQ, 1280, ‘Window’, window, ‘OverlapLength’, overlap_length, ‘FFTLength’, fft_points);
new_data_Set_STFT{index_onearray,2} = F; %Frequency
% Amplitude of Spectogram
new_data_Set_STFT{index_onearray,3} =mean(abs(S’));
% IQ data phase
new_data_Set_STFT{index_onearray,4} = atan(imag(S)/real(S));
FFT:
% IQ DATA FFT
new_data_Set{index_onearray,2} = fft(IQ(:,index_global));
% Frequency length
N = length(IQ(:,index_global)); % FFT length
fs = 256; % sampling frequency
calculation = (0:N-1)*(fs/N);
f1 = reshape(calculation, 1280, …
1);
new_data_Set{index_onearray,3} = f1; % Frequency vector
% Amplitude of Spectogram
new_data_Set{index_onearray,4} = abs(IQ(:,index_global));
% IQ data phase
new_data_Set{index_onearray,5} = angle(IQ(:,index_global));
After Plotting, STFT spektogram looks like this. In frequency and time(x time, y frequency) scale it looks wrong. What’s wrong with my code?
FFT:
Thanks for help. radar, signal processing, spectogram, fft, stft, fmcw, plotting, iq data MATLAB Answers — New Questions
solve using trapezoidal rule
Write a local function named integration to numerically integrate the tabulated values by using trapezoidal rule, if X and Y values are user input and difference between any two consecutive X values is same. Use the above local function to determine ∫ 𝑦 𝑑𝑥 from 0 to π by using following data. X= 0, π/6 ,π/3 ,π/2, 2π/3, 5π/6, π and Y= 0, 0.5 ,0.866, 1 ,0.866, 0.5, 0 .Validate the result by using appropriate built-in functionWrite a local function named integration to numerically integrate the tabulated values by using trapezoidal rule, if X and Y values are user input and difference between any two consecutive X values is same. Use the above local function to determine ∫ 𝑦 𝑑𝑥 from 0 to π by using following data. X= 0, π/6 ,π/3 ,π/2, 2π/3, 5π/6, π and Y= 0, 0.5 ,0.866, 1 ,0.866, 0.5, 0 .Validate the result by using appropriate built-in function Write a local function named integration to numerically integrate the tabulated values by using trapezoidal rule, if X and Y values are user input and difference between any two consecutive X values is same. Use the above local function to determine ∫ 𝑦 𝑑𝑥 from 0 to π by using following data. X= 0, π/6 ,π/3 ,π/2, 2π/3, 5π/6, π and Y= 0, 0.5 ,0.866, 1 ,0.866, 0.5, 0 .Validate the result by using appropriate built-in function optimization, trapezoidal rule MATLAB Answers — New Questions
error when using {}
Hi guys, im trying to run a software by Matlab and do some analysis for some models and this ‘Periodandfrequencies’ is an example of some results. But I am facing error when I want to keep the result inside the cell-array. Basicaly I want it to be iix1 cell aray and each cell is 12×7 Table.Hi guys, im trying to run a software by Matlab and do some analysis for some models and this ‘Periodandfrequencies’ is an example of some results. But I am facing error when I want to keep the result inside the cell-array. Basicaly I want it to be iix1 cell aray and each cell is 12×7 Table. Hi guys, im trying to run a software by Matlab and do some analysis for some models and this ‘Periodandfrequencies’ is an example of some results. But I am facing error when I want to keep the result inside the cell-array. Basicaly I want it to be iix1 cell aray and each cell is 12×7 Table. cell array MATLAB Answers — New Questions
How to evaluate Trained detector on test dataset?
Hello,
How can I evaluate Trained detector on test dataset, pelase? Here example code:
tic
detectionResults = detect(Trained_detector, testData, ‘MiniBatchSize’, 50);
% Evaluate the object detector using average precision metric.
[ap,recall,precision] = evaluateDetectionPrecision(detectionResults, testData);
classID = 1;
figure
plot(recall{classID},precision{classID},’LineWidth’,4)
xlabel("recall")
ylabel("accuracy")
grid on
title(sprintf("AP = %.3f",ap(classID)))
toc
Here Error:
Expected I to be one of these types:
double, single, int16, uint16, uint8, logical
Instead its type was matlab.io.datastore.CombinedDatastore.Hello,
How can I evaluate Trained detector on test dataset, pelase? Here example code:
tic
detectionResults = detect(Trained_detector, testData, ‘MiniBatchSize’, 50);
% Evaluate the object detector using average precision metric.
[ap,recall,precision] = evaluateDetectionPrecision(detectionResults, testData);
classID = 1;
figure
plot(recall{classID},precision{classID},’LineWidth’,4)
xlabel("recall")
ylabel("accuracy")
grid on
title(sprintf("AP = %.3f",ap(classID)))
toc
Here Error:
Expected I to be one of these types:
double, single, int16, uint16, uint8, logical
Instead its type was matlab.io.datastore.CombinedDatastore. Hello,
How can I evaluate Trained detector on test dataset, pelase? Here example code:
tic
detectionResults = detect(Trained_detector, testData, ‘MiniBatchSize’, 50);
% Evaluate the object detector using average precision metric.
[ap,recall,precision] = evaluateDetectionPrecision(detectionResults, testData);
classID = 1;
figure
plot(recall{classID},precision{classID},’LineWidth’,4)
xlabel("recall")
ylabel("accuracy")
grid on
title(sprintf("AP = %.3f",ap(classID)))
toc
Here Error:
Expected I to be one of these types:
double, single, int16, uint16, uint8, logical
Instead its type was matlab.io.datastore.CombinedDatastore. image processing, deep learning, image analysis, computer vision MATLAB Answers — New Questions
Need help solving heat equation using adi method
This is my equation, which I solve using the ADI scheme, but it doesn’t work correctly; it only works with a very small dt. Can anyone help me?
Nx = 50;
Ny = 50;
dx = 1.0 / (Nx – 1);
dy = 2.0 / (Ny – 1);
T = 0.00001;
dt = 0.000000001;
Nt = 1000; % Total number of time steps
alpha = 4; % Diffusion coefficient
x = linspace(0, 1, Nx);
y = linspace(0, 2, Ny);
[X, Y] = meshgrid(x, y);
U = Y .* X + 1; % Initial condition
function val = f(x, y, t)
val = exp(t) * cos(pi * x / 2) * sin(pi * y / 4);
end
function U = apply_boundary_conditions(U, x, y, dy, dx)
U(:, 1) = 1; % y=0
U(:, end) = U(:, end-1) + dy .* x’; % y=2
U(1, 🙂 = U(2, 🙂 – dx * y; % x=0
U(end, 🙂 = y’ + 1; % x=1
end
% Thomas algorithm
function x = thomas_algorithm(a, b, c, d)
n = length(d);
c_star = zeros(n-1, 1);
d_star = zeros(n, 1);
x = zeros(n, 1);
c_star(1) = c(1) / b(1);
d_star(1) = d(1) / b(1);
for i = 2:n-1
temp = b(i) – a(i-1) * c_star(i-1);
c_star(i) = c(i) / temp;
d_star(i) = (d(i) – a(i-1) * d_star(i-1)) / temp;
end
d_star(n) = (d(n) – a(n-1) * d_star(n-1)) / b(n);
x(n) = d_star(n);
for i = n-1:-1:1
x(i) = d_star(i) – c_star(i) * x(i+1);
end
end
% ADI method
for n = 1:Nt
U = apply_boundary_conditions(U, x, y, dy, dx);
% First half-step: X-direction implicit, Y-direction explicit
for j = 2:Ny-1
a = -alpha * dt / (2 * dx^2) * ones(Nx-1, 1);
b = (1 + alpha * dt / dx^2) * ones(Nx, 1);
c = -alpha * dt / (2 * dx^2) * ones(Nx-1, 1);
d = U(j, 2:end-1)’ + 0.5 * alpha * dt / dy^2 * (U(j+1, 2:end-1) – 2 * U(j, 2:end-1) + U(j-1, 2:end-1))’ + dt * f(x(2:end-1), y(j), n*dt);
U(j, 2:end-1) = thomas_algorithm(a, b(2:end-1), c, d);
end
% Second half-step: Y-direction implicit, X-direction explicit
for i = 2:Nx-1
a = -alpha * dt / (2 * dy^2) * ones(Ny-1, 1);
b = (1 + alpha * dt / dy^2) * ones(Ny, 1);
c = -alpha * dt / (2 * dy^2) * ones(Ny-1, 1);
d = U(2:end-1, i) + 0.5 * alpha * dt / dx^2 * (U(2:end-1, i+1) – 2 * U(2:end-1, i) + U(2:end-1, i-1)) + dt * f(x(i), y(2:end-1), n*dt);
U(2:end-1, i) = thomas_algorithm(a, b(2:end-1), c, d);
end
U = apply_boundary_conditions(U, x, y, dy, dx);
end
% Visualization
surf(X, Y, U);
title(‘ADI Solution at T=’ + string(T));
xlabel(‘X’);
ylabel(‘Y’);
zlabel(‘U’);
colorbar;This is my equation, which I solve using the ADI scheme, but it doesn’t work correctly; it only works with a very small dt. Can anyone help me?
Nx = 50;
Ny = 50;
dx = 1.0 / (Nx – 1);
dy = 2.0 / (Ny – 1);
T = 0.00001;
dt = 0.000000001;
Nt = 1000; % Total number of time steps
alpha = 4; % Diffusion coefficient
x = linspace(0, 1, Nx);
y = linspace(0, 2, Ny);
[X, Y] = meshgrid(x, y);
U = Y .* X + 1; % Initial condition
function val = f(x, y, t)
val = exp(t) * cos(pi * x / 2) * sin(pi * y / 4);
end
function U = apply_boundary_conditions(U, x, y, dy, dx)
U(:, 1) = 1; % y=0
U(:, end) = U(:, end-1) + dy .* x’; % y=2
U(1, 🙂 = U(2, 🙂 – dx * y; % x=0
U(end, 🙂 = y’ + 1; % x=1
end
% Thomas algorithm
function x = thomas_algorithm(a, b, c, d)
n = length(d);
c_star = zeros(n-1, 1);
d_star = zeros(n, 1);
x = zeros(n, 1);
c_star(1) = c(1) / b(1);
d_star(1) = d(1) / b(1);
for i = 2:n-1
temp = b(i) – a(i-1) * c_star(i-1);
c_star(i) = c(i) / temp;
d_star(i) = (d(i) – a(i-1) * d_star(i-1)) / temp;
end
d_star(n) = (d(n) – a(n-1) * d_star(n-1)) / b(n);
x(n) = d_star(n);
for i = n-1:-1:1
x(i) = d_star(i) – c_star(i) * x(i+1);
end
end
% ADI method
for n = 1:Nt
U = apply_boundary_conditions(U, x, y, dy, dx);
% First half-step: X-direction implicit, Y-direction explicit
for j = 2:Ny-1
a = -alpha * dt / (2 * dx^2) * ones(Nx-1, 1);
b = (1 + alpha * dt / dx^2) * ones(Nx, 1);
c = -alpha * dt / (2 * dx^2) * ones(Nx-1, 1);
d = U(j, 2:end-1)’ + 0.5 * alpha * dt / dy^2 * (U(j+1, 2:end-1) – 2 * U(j, 2:end-1) + U(j-1, 2:end-1))’ + dt * f(x(2:end-1), y(j), n*dt);
U(j, 2:end-1) = thomas_algorithm(a, b(2:end-1), c, d);
end
% Second half-step: Y-direction implicit, X-direction explicit
for i = 2:Nx-1
a = -alpha * dt / (2 * dy^2) * ones(Ny-1, 1);
b = (1 + alpha * dt / dy^2) * ones(Ny, 1);
c = -alpha * dt / (2 * dy^2) * ones(Ny-1, 1);
d = U(2:end-1, i) + 0.5 * alpha * dt / dx^2 * (U(2:end-1, i+1) – 2 * U(2:end-1, i) + U(2:end-1, i-1)) + dt * f(x(i), y(2:end-1), n*dt);
U(2:end-1, i) = thomas_algorithm(a, b(2:end-1), c, d);
end
U = apply_boundary_conditions(U, x, y, dy, dx);
end
% Visualization
surf(X, Y, U);
title(‘ADI Solution at T=’ + string(T));
xlabel(‘X’);
ylabel(‘Y’);
zlabel(‘U’);
colorbar; This is my equation, which I solve using the ADI scheme, but it doesn’t work correctly; it only works with a very small dt. Can anyone help me?
Nx = 50;
Ny = 50;
dx = 1.0 / (Nx – 1);
dy = 2.0 / (Ny – 1);
T = 0.00001;
dt = 0.000000001;
Nt = 1000; % Total number of time steps
alpha = 4; % Diffusion coefficient
x = linspace(0, 1, Nx);
y = linspace(0, 2, Ny);
[X, Y] = meshgrid(x, y);
U = Y .* X + 1; % Initial condition
function val = f(x, y, t)
val = exp(t) * cos(pi * x / 2) * sin(pi * y / 4);
end
function U = apply_boundary_conditions(U, x, y, dy, dx)
U(:, 1) = 1; % y=0
U(:, end) = U(:, end-1) + dy .* x’; % y=2
U(1, 🙂 = U(2, 🙂 – dx * y; % x=0
U(end, 🙂 = y’ + 1; % x=1
end
% Thomas algorithm
function x = thomas_algorithm(a, b, c, d)
n = length(d);
c_star = zeros(n-1, 1);
d_star = zeros(n, 1);
x = zeros(n, 1);
c_star(1) = c(1) / b(1);
d_star(1) = d(1) / b(1);
for i = 2:n-1
temp = b(i) – a(i-1) * c_star(i-1);
c_star(i) = c(i) / temp;
d_star(i) = (d(i) – a(i-1) * d_star(i-1)) / temp;
end
d_star(n) = (d(n) – a(n-1) * d_star(n-1)) / b(n);
x(n) = d_star(n);
for i = n-1:-1:1
x(i) = d_star(i) – c_star(i) * x(i+1);
end
end
% ADI method
for n = 1:Nt
U = apply_boundary_conditions(U, x, y, dy, dx);
% First half-step: X-direction implicit, Y-direction explicit
for j = 2:Ny-1
a = -alpha * dt / (2 * dx^2) * ones(Nx-1, 1);
b = (1 + alpha * dt / dx^2) * ones(Nx, 1);
c = -alpha * dt / (2 * dx^2) * ones(Nx-1, 1);
d = U(j, 2:end-1)’ + 0.5 * alpha * dt / dy^2 * (U(j+1, 2:end-1) – 2 * U(j, 2:end-1) + U(j-1, 2:end-1))’ + dt * f(x(2:end-1), y(j), n*dt);
U(j, 2:end-1) = thomas_algorithm(a, b(2:end-1), c, d);
end
% Second half-step: Y-direction implicit, X-direction explicit
for i = 2:Nx-1
a = -alpha * dt / (2 * dy^2) * ones(Ny-1, 1);
b = (1 + alpha * dt / dy^2) * ones(Ny, 1);
c = -alpha * dt / (2 * dy^2) * ones(Ny-1, 1);
d = U(2:end-1, i) + 0.5 * alpha * dt / dx^2 * (U(2:end-1, i+1) – 2 * U(2:end-1, i) + U(2:end-1, i-1)) + dt * f(x(i), y(2:end-1), n*dt);
U(2:end-1, i) = thomas_algorithm(a, b(2:end-1), c, d);
end
U = apply_boundary_conditions(U, x, y, dy, dx);
end
% Visualization
surf(X, Y, U);
title(‘ADI Solution at T=’ + string(T));
xlabel(‘X’);
ylabel(‘Y’);
zlabel(‘U’);
colorbar; adischeme heat equations MATLAB Answers — New Questions
Parfor “Out of Memory during deserialization” in large matrix operations
I am trying to use parfor for a large matrix operation. I am getting Out of Memory during deserialization error. Is there a way to minimize the memory? Below is the minimal example code:
clc; clear;
warpedImages = num2cell(uint8(randi([0,255], 1654, 6288, 3, 35)),1:3);
warpedImages = reshape(warpedImages,1,[]);
% Initialze
n = length(warpedImages);
sigmaN = 10;
sigmag = 0.1;
panoramasize = size(warpedImages{1});
Amat = cell(n);
Bvec = zeros(n,1);
IuppeIdx = nonzeros(triu(reshape(1:numel(Amat), size(Amat))));
Amat_temp = cell(1,length(IuppeIdx));
matSize = size(Amat);
% 4D warped images (Slicing)
wim_4d = cell2mat(reshape(warpedImages,1,1,1,[]));
% Get the Ibarijs and Nijs
parfor i = 1:length(IuppeIdx)
% Index to subscripts
[ii,jj] = ind2sub(matSize, IuppeIdx(i));
if ii == jj
diag_val_1 = 0;
diag_val_2 = 0;
Z = 1:n;
Z(Z==ii) = [];
for d = Z
[Ibarij, Ibarji, Nij] = getIbarNij(panoramasize, wim_4d(:,:,:,ii), wim_4d(:,:,:,d));
diag_val_1 = diag_val_1 + ( (Nij + Nij) .* Ibarij.^2 );
diag_val_2 = diag_val_2 + Nij;
end
diag_val = diag_val_1 + (sigmaN^2/sigmag^2) * diag_val_2;
B_val = (sigmaN^2/sigmag^2) * diag_val_2;
Amat_temp{i} = diag_val;
Bvec(i) = B_val
end
if ii ~= jj
[Ibarij,Ibarji,Nij] = getIbarNij(panoramasize, wim_4d(:,:,:,ii), wim_4d(:,:,:,jj));
Amat_temp{i} = -(Nij+Nij) .* (Ibarij .* Ibarji);
end
end
function [Ibarij,Ibarji,Nij] = getIbarNij(panoramasize, Imij, Imji)
Ibarij = zeros(panoramasize,’uint8′);
Ibarji = zeros(panoramasize,’uint8′);
% Overlay the warpedImage onto the panorama.
maski = imbinarize(rgb2gray(255 * Imij));
maskj = imbinarize(rgb2gray(255 * Imji));
% Find the overlap mask
Nij_im = maski & maskj;
Nij_im = imfill(Nij_im, ‘holes’);
Nijidx = repmat(Nij_im, 1, 1, size(Imij,3));
% Get the overlapping region RGB values for two images
Ibarij(Nijidx) = Imij(Nijidx);
Ibarji(Nijidx) = Imji(Nijidx);
% Convert to double
Ibarij_double = double(Ibarij);
Ibarji_double = double(Ibarji);
% Nij
Nij = sum(sum(Nij_im));
% Ibar ijs
Ibarij = reshape(sum(sum(Ibarij_double)) ./ Nij, 1, 3);
Ibarji = reshape(sum(sum(Ibarji_double)) ./ Nij, 1, 3);
% Replace NaNs by zeros
Ibarij(isnan(Ibarij)) = 0;
Ibarji(isnan(Ibarji)) = 0;
end
Line [Ibarij, Ibarji, Nij] = getIbarNij(panoramasize, wim_4d(:,:,:,ii), wim_4d(:,:,:,d)); throws a warning message: The entire array or structure ‘wim_4d’ is a broadcast variable. This might result in unnecessary communication overhead. I have used ind2sub for getting the subscripts as it is easy to work. However, wim_4d(:,:,:,ii) and others cannot be sliced. Any other suggestion and help is appreciated!I am trying to use parfor for a large matrix operation. I am getting Out of Memory during deserialization error. Is there a way to minimize the memory? Below is the minimal example code:
clc; clear;
warpedImages = num2cell(uint8(randi([0,255], 1654, 6288, 3, 35)),1:3);
warpedImages = reshape(warpedImages,1,[]);
% Initialze
n = length(warpedImages);
sigmaN = 10;
sigmag = 0.1;
panoramasize = size(warpedImages{1});
Amat = cell(n);
Bvec = zeros(n,1);
IuppeIdx = nonzeros(triu(reshape(1:numel(Amat), size(Amat))));
Amat_temp = cell(1,length(IuppeIdx));
matSize = size(Amat);
% 4D warped images (Slicing)
wim_4d = cell2mat(reshape(warpedImages,1,1,1,[]));
% Get the Ibarijs and Nijs
parfor i = 1:length(IuppeIdx)
% Index to subscripts
[ii,jj] = ind2sub(matSize, IuppeIdx(i));
if ii == jj
diag_val_1 = 0;
diag_val_2 = 0;
Z = 1:n;
Z(Z==ii) = [];
for d = Z
[Ibarij, Ibarji, Nij] = getIbarNij(panoramasize, wim_4d(:,:,:,ii), wim_4d(:,:,:,d));
diag_val_1 = diag_val_1 + ( (Nij + Nij) .* Ibarij.^2 );
diag_val_2 = diag_val_2 + Nij;
end
diag_val = diag_val_1 + (sigmaN^2/sigmag^2) * diag_val_2;
B_val = (sigmaN^2/sigmag^2) * diag_val_2;
Amat_temp{i} = diag_val;
Bvec(i) = B_val
end
if ii ~= jj
[Ibarij,Ibarji,Nij] = getIbarNij(panoramasize, wim_4d(:,:,:,ii), wim_4d(:,:,:,jj));
Amat_temp{i} = -(Nij+Nij) .* (Ibarij .* Ibarji);
end
end
function [Ibarij,Ibarji,Nij] = getIbarNij(panoramasize, Imij, Imji)
Ibarij = zeros(panoramasize,’uint8′);
Ibarji = zeros(panoramasize,’uint8′);
% Overlay the warpedImage onto the panorama.
maski = imbinarize(rgb2gray(255 * Imij));
maskj = imbinarize(rgb2gray(255 * Imji));
% Find the overlap mask
Nij_im = maski & maskj;
Nij_im = imfill(Nij_im, ‘holes’);
Nijidx = repmat(Nij_im, 1, 1, size(Imij,3));
% Get the overlapping region RGB values for two images
Ibarij(Nijidx) = Imij(Nijidx);
Ibarji(Nijidx) = Imji(Nijidx);
% Convert to double
Ibarij_double = double(Ibarij);
Ibarji_double = double(Ibarji);
% Nij
Nij = sum(sum(Nij_im));
% Ibar ijs
Ibarij = reshape(sum(sum(Ibarij_double)) ./ Nij, 1, 3);
Ibarji = reshape(sum(sum(Ibarji_double)) ./ Nij, 1, 3);
% Replace NaNs by zeros
Ibarij(isnan(Ibarij)) = 0;
Ibarji(isnan(Ibarji)) = 0;
end
Line [Ibarij, Ibarji, Nij] = getIbarNij(panoramasize, wim_4d(:,:,:,ii), wim_4d(:,:,:,d)); throws a warning message: The entire array or structure ‘wim_4d’ is a broadcast variable. This might result in unnecessary communication overhead. I have used ind2sub for getting the subscripts as it is easy to work. However, wim_4d(:,:,:,ii) and others cannot be sliced. Any other suggestion and help is appreciated! I am trying to use parfor for a large matrix operation. I am getting Out of Memory during deserialization error. Is there a way to minimize the memory? Below is the minimal example code:
clc; clear;
warpedImages = num2cell(uint8(randi([0,255], 1654, 6288, 3, 35)),1:3);
warpedImages = reshape(warpedImages,1,[]);
% Initialze
n = length(warpedImages);
sigmaN = 10;
sigmag = 0.1;
panoramasize = size(warpedImages{1});
Amat = cell(n);
Bvec = zeros(n,1);
IuppeIdx = nonzeros(triu(reshape(1:numel(Amat), size(Amat))));
Amat_temp = cell(1,length(IuppeIdx));
matSize = size(Amat);
% 4D warped images (Slicing)
wim_4d = cell2mat(reshape(warpedImages,1,1,1,[]));
% Get the Ibarijs and Nijs
parfor i = 1:length(IuppeIdx)
% Index to subscripts
[ii,jj] = ind2sub(matSize, IuppeIdx(i));
if ii == jj
diag_val_1 = 0;
diag_val_2 = 0;
Z = 1:n;
Z(Z==ii) = [];
for d = Z
[Ibarij, Ibarji, Nij] = getIbarNij(panoramasize, wim_4d(:,:,:,ii), wim_4d(:,:,:,d));
diag_val_1 = diag_val_1 + ( (Nij + Nij) .* Ibarij.^2 );
diag_val_2 = diag_val_2 + Nij;
end
diag_val = diag_val_1 + (sigmaN^2/sigmag^2) * diag_val_2;
B_val = (sigmaN^2/sigmag^2) * diag_val_2;
Amat_temp{i} = diag_val;
Bvec(i) = B_val
end
if ii ~= jj
[Ibarij,Ibarji,Nij] = getIbarNij(panoramasize, wim_4d(:,:,:,ii), wim_4d(:,:,:,jj));
Amat_temp{i} = -(Nij+Nij) .* (Ibarij .* Ibarji);
end
end
function [Ibarij,Ibarji,Nij] = getIbarNij(panoramasize, Imij, Imji)
Ibarij = zeros(panoramasize,’uint8′);
Ibarji = zeros(panoramasize,’uint8′);
% Overlay the warpedImage onto the panorama.
maski = imbinarize(rgb2gray(255 * Imij));
maskj = imbinarize(rgb2gray(255 * Imji));
% Find the overlap mask
Nij_im = maski & maskj;
Nij_im = imfill(Nij_im, ‘holes’);
Nijidx = repmat(Nij_im, 1, 1, size(Imij,3));
% Get the overlapping region RGB values for two images
Ibarij(Nijidx) = Imij(Nijidx);
Ibarji(Nijidx) = Imji(Nijidx);
% Convert to double
Ibarij_double = double(Ibarij);
Ibarji_double = double(Ibarji);
% Nij
Nij = sum(sum(Nij_im));
% Ibar ijs
Ibarij = reshape(sum(sum(Ibarij_double)) ./ Nij, 1, 3);
Ibarji = reshape(sum(sum(Ibarji_double)) ./ Nij, 1, 3);
% Replace NaNs by zeros
Ibarij(isnan(Ibarij)) = 0;
Ibarji(isnan(Ibarji)) = 0;
end
Line [Ibarij, Ibarji, Nij] = getIbarNij(panoramasize, wim_4d(:,:,:,ii), wim_4d(:,:,:,d)); throws a warning message: The entire array or structure ‘wim_4d’ is a broadcast variable. This might result in unnecessary communication overhead. I have used ind2sub for getting the subscripts as it is easy to work. However, wim_4d(:,:,:,ii) and others cannot be sliced. Any other suggestion and help is appreciated! parfor, memory MATLAB Answers — New Questions
How do I get my old university account as in gmail account?
I have one old account with university email address but now I dont have that email address so I need to connect that account to my new gmail account so How can I do it?
I dont have access to my old account anymore. Thank you
old account
https://in.mathworks.com/matlabcentral/profile/authors/14080250
New account is
https://in.mathworks.com/matlabcentral/profile/authors/14080250I have one old account with university email address but now I dont have that email address so I need to connect that account to my new gmail account so How can I do it?
I dont have access to my old account anymore. Thank you
old account
https://in.mathworks.com/matlabcentral/profile/authors/14080250
New account is
https://in.mathworks.com/matlabcentral/profile/authors/14080250 I have one old account with university email address but now I dont have that email address so I need to connect that account to my new gmail account so How can I do it?
I dont have access to my old account anymore. Thank you
old account
https://in.mathworks.com/matlabcentral/profile/authors/14080250
New account is
https://in.mathworks.com/matlabcentral/profile/authors/14080250 support, matlab, account MATLAB Answers — New Questions
I need help with Predictor-Corrector Scheme for PDE. I use 2d matrices instead of 3d matrices, but something is wrong and I don’t understand what.
This is my problem.
∂u/∂t = 4∆u + e^t * cos(πx/2) * sin(πy/4) x ∈ (0, 1), y ∈ (0, 2), t ∈ (0, T],
∂u/∂x | x=0 = y,
u | x=1 = y + 1,
u | y=0 = 1,
∂u/∂y | y=2 = x,
u | t=0 = x*y + 1
% CODE
% Parameters
Lx = 1.0; Ly = 2.0;
Nx = 50; Ny = 100; % Number of steps in x and y directions
% Coefficient
D = 4;
% Time parameters
T = 0.0001;
Nt = 1000; % Number of time steps
dt = T / Nt;
% Grid
dx = Lx / (Nx – 0.5);
x = linspace(-dx / 2, Lx, Nx + 1);
dy = Ly / (Ny – 0.5);
y = linspace(Ly + dy / 2, 0, Ny + 1);
t = linspace(0, T, Nt + 1);
% Initial condition
[X, Y] = meshgrid(x, y);
u = X .* Y + 1;
% Inhomogeneity function
f = @(x, y, t) exp(t) .* cos(pi * x / 2) .* sin(pi * y / 4);
% Thomas algorithm for tridiagonal system
function x = thomas_algorithm(a, b, c, d)
n = length(d);
cp = zeros(n – 1, 1);
dp = zeros(n, 1);
% Forward sweep
cp(1) = c(1) / b(1);
dp(1) = d(1) / b(1);
for i = 2:n
m = 1.0 / (b(i) – a(i – 1) * cp(i – 1));
if i < n
cp(i) = c(i) * m;
end
dp(i) = (d(i) – a(i – 1) * dp(i – 1)) * m;
end
% Backward substitution
x = zeros(n, 1);
x(n) = dp(n);
for i = n – 1:-1:1
x(i) = dp(i) – cp(i) * x(i + 1);
end
end
% Numerical solution scheme
for n = 1:Nt
u_old = u;
% Predictor step in x direction
for j = 2:Ny
a = -D * ones(Nx – 1, 1) / (dx^2);
b = (1 / (dt / 2) + 2 * D / (dx^2)) * ones(Nx – 1, 1);
c = -D * ones(Nx – 1, 1) / (dx^2);
d = u(2:Nx, j) / (dt / 2) + f(x(2:Nx), y(j), t(n));
u(2:Nx, j) = thomas_algorithm(a, b, c, d);
end
% Predictor step in y direction
for i = 2:Nx
a = -D * ones(Ny – 1, 1) / (dy^2);
b = (1 / (dt / 2) + 2 * D / (dy^2)) * ones(Ny – 1, 1);
c = -D * ones(Ny – 1, 1) / (dy^2);
d = u(i, 2:Ny) / (dt / 2) + f(x(i), y(2:Ny), t(n));
u(i, 2:Ny) = thomas_algorithm(a, b, c, d);
end
u_old1 = u;
% Corrector step
for i = 2:Nx
for j = 2:Ny
u(i, j) = u_old(i, j) + dt * (…
D * (u_old1(i – 1, j) – 2 * u_old1(i, j) + u_old1(i + 1, j)) / (dx^2) + …
D * (u_old1(i, j – 1) – 2 * u_old1(i, j) + u_old1(i, j + 1)) / (dy^2) + …
f(x(i), y(j), t(n) + dt / 2) …
);
end
end
% Apply boundary conditions
u(:, 1) = 1; % y = 0
u(:, end) = u(:, end – 1) + dy * x; % Neumann y = 2
u(1, 🙂 = u(2, 🙂 – dx * y; % Neumann x = 0
u(end, 🙂 = y + 1; % x = 1
endThis is my problem.
∂u/∂t = 4∆u + e^t * cos(πx/2) * sin(πy/4) x ∈ (0, 1), y ∈ (0, 2), t ∈ (0, T],
∂u/∂x | x=0 = y,
u | x=1 = y + 1,
u | y=0 = 1,
∂u/∂y | y=2 = x,
u | t=0 = x*y + 1
% CODE
% Parameters
Lx = 1.0; Ly = 2.0;
Nx = 50; Ny = 100; % Number of steps in x and y directions
% Coefficient
D = 4;
% Time parameters
T = 0.0001;
Nt = 1000; % Number of time steps
dt = T / Nt;
% Grid
dx = Lx / (Nx – 0.5);
x = linspace(-dx / 2, Lx, Nx + 1);
dy = Ly / (Ny – 0.5);
y = linspace(Ly + dy / 2, 0, Ny + 1);
t = linspace(0, T, Nt + 1);
% Initial condition
[X, Y] = meshgrid(x, y);
u = X .* Y + 1;
% Inhomogeneity function
f = @(x, y, t) exp(t) .* cos(pi * x / 2) .* sin(pi * y / 4);
% Thomas algorithm for tridiagonal system
function x = thomas_algorithm(a, b, c, d)
n = length(d);
cp = zeros(n – 1, 1);
dp = zeros(n, 1);
% Forward sweep
cp(1) = c(1) / b(1);
dp(1) = d(1) / b(1);
for i = 2:n
m = 1.0 / (b(i) – a(i – 1) * cp(i – 1));
if i < n
cp(i) = c(i) * m;
end
dp(i) = (d(i) – a(i – 1) * dp(i – 1)) * m;
end
% Backward substitution
x = zeros(n, 1);
x(n) = dp(n);
for i = n – 1:-1:1
x(i) = dp(i) – cp(i) * x(i + 1);
end
end
% Numerical solution scheme
for n = 1:Nt
u_old = u;
% Predictor step in x direction
for j = 2:Ny
a = -D * ones(Nx – 1, 1) / (dx^2);
b = (1 / (dt / 2) + 2 * D / (dx^2)) * ones(Nx – 1, 1);
c = -D * ones(Nx – 1, 1) / (dx^2);
d = u(2:Nx, j) / (dt / 2) + f(x(2:Nx), y(j), t(n));
u(2:Nx, j) = thomas_algorithm(a, b, c, d);
end
% Predictor step in y direction
for i = 2:Nx
a = -D * ones(Ny – 1, 1) / (dy^2);
b = (1 / (dt / 2) + 2 * D / (dy^2)) * ones(Ny – 1, 1);
c = -D * ones(Ny – 1, 1) / (dy^2);
d = u(i, 2:Ny) / (dt / 2) + f(x(i), y(2:Ny), t(n));
u(i, 2:Ny) = thomas_algorithm(a, b, c, d);
end
u_old1 = u;
% Corrector step
for i = 2:Nx
for j = 2:Ny
u(i, j) = u_old(i, j) + dt * (…
D * (u_old1(i – 1, j) – 2 * u_old1(i, j) + u_old1(i + 1, j)) / (dx^2) + …
D * (u_old1(i, j – 1) – 2 * u_old1(i, j) + u_old1(i, j + 1)) / (dy^2) + …
f(x(i), y(j), t(n) + dt / 2) …
);
end
end
% Apply boundary conditions
u(:, 1) = 1; % y = 0
u(:, end) = u(:, end – 1) + dy * x; % Neumann y = 2
u(1, 🙂 = u(2, 🙂 – dx * y; % Neumann x = 0
u(end, 🙂 = y + 1; % x = 1
end This is my problem.
∂u/∂t = 4∆u + e^t * cos(πx/2) * sin(πy/4) x ∈ (0, 1), y ∈ (0, 2), t ∈ (0, T],
∂u/∂x | x=0 = y,
u | x=1 = y + 1,
u | y=0 = 1,
∂u/∂y | y=2 = x,
u | t=0 = x*y + 1
% CODE
% Parameters
Lx = 1.0; Ly = 2.0;
Nx = 50; Ny = 100; % Number of steps in x and y directions
% Coefficient
D = 4;
% Time parameters
T = 0.0001;
Nt = 1000; % Number of time steps
dt = T / Nt;
% Grid
dx = Lx / (Nx – 0.5);
x = linspace(-dx / 2, Lx, Nx + 1);
dy = Ly / (Ny – 0.5);
y = linspace(Ly + dy / 2, 0, Ny + 1);
t = linspace(0, T, Nt + 1);
% Initial condition
[X, Y] = meshgrid(x, y);
u = X .* Y + 1;
% Inhomogeneity function
f = @(x, y, t) exp(t) .* cos(pi * x / 2) .* sin(pi * y / 4);
% Thomas algorithm for tridiagonal system
function x = thomas_algorithm(a, b, c, d)
n = length(d);
cp = zeros(n – 1, 1);
dp = zeros(n, 1);
% Forward sweep
cp(1) = c(1) / b(1);
dp(1) = d(1) / b(1);
for i = 2:n
m = 1.0 / (b(i) – a(i – 1) * cp(i – 1));
if i < n
cp(i) = c(i) * m;
end
dp(i) = (d(i) – a(i – 1) * dp(i – 1)) * m;
end
% Backward substitution
x = zeros(n, 1);
x(n) = dp(n);
for i = n – 1:-1:1
x(i) = dp(i) – cp(i) * x(i + 1);
end
end
% Numerical solution scheme
for n = 1:Nt
u_old = u;
% Predictor step in x direction
for j = 2:Ny
a = -D * ones(Nx – 1, 1) / (dx^2);
b = (1 / (dt / 2) + 2 * D / (dx^2)) * ones(Nx – 1, 1);
c = -D * ones(Nx – 1, 1) / (dx^2);
d = u(2:Nx, j) / (dt / 2) + f(x(2:Nx), y(j), t(n));
u(2:Nx, j) = thomas_algorithm(a, b, c, d);
end
% Predictor step in y direction
for i = 2:Nx
a = -D * ones(Ny – 1, 1) / (dy^2);
b = (1 / (dt / 2) + 2 * D / (dy^2)) * ones(Ny – 1, 1);
c = -D * ones(Ny – 1, 1) / (dy^2);
d = u(i, 2:Ny) / (dt / 2) + f(x(i), y(2:Ny), t(n));
u(i, 2:Ny) = thomas_algorithm(a, b, c, d);
end
u_old1 = u;
% Corrector step
for i = 2:Nx
for j = 2:Ny
u(i, j) = u_old(i, j) + dt * (…
D * (u_old1(i – 1, j) – 2 * u_old1(i, j) + u_old1(i + 1, j)) / (dx^2) + …
D * (u_old1(i, j – 1) – 2 * u_old1(i, j) + u_old1(i, j + 1)) / (dy^2) + …
f(x(i), y(j), t(n) + dt / 2) …
);
end
end
% Apply boundary conditions
u(:, 1) = 1; % y = 0
u(:, end) = u(:, end – 1) + dy * x; % Neumann y = 2
u(1, 🙂 = u(2, 🙂 – dx * y; % Neumann x = 0
u(end, 🙂 = y + 1; % x = 1
end pde, predictor-corrector, 2d MATLAB Answers — New Questions
Stateflow output string error
"String buffer overflow. Increase the maximum length of the output string data type."
this is the error showing up as soon as try to run my simulink model which contains stateflow chart.
the stateflow chart contains some strings that needs to be the output and can anybody tell me how can I increase the max length."String buffer overflow. Increase the maximum length of the output string data type."
this is the error showing up as soon as try to run my simulink model which contains stateflow chart.
the stateflow chart contains some strings that needs to be the output and can anybody tell me how can I increase the max length. "String buffer overflow. Increase the maximum length of the output string data type."
this is the error showing up as soon as try to run my simulink model which contains stateflow chart.
the stateflow chart contains some strings that needs to be the output and can anybody tell me how can I increase the max length. simulink, stateflow, matlab, string MATLAB Answers — New Questions
how to run MLP
HI.EVER
my teacher give me tasks and i’m newbi in this field…i want to estimaition this function by (k_means and rbf) method
first step: i need to learn how to write function like it on matlab
seconed step : give me the help i can understand how to run mlp or rbfHI.EVER
my teacher give me tasks and i’m newbi in this field…i want to estimaition this function by (k_means and rbf) method
first step: i need to learn how to write function like it on matlab
seconed step : give me the help i can understand how to run mlp or rbf HI.EVER
my teacher give me tasks and i’m newbi in this field…i want to estimaition this function by (k_means and rbf) method
first step: i need to learn how to write function like it on matlab
seconed step : give me the help i can understand how to run mlp or rbf mlp, rbf MATLAB Answers — New Questions
Matlab Home. Possible to get the coders?
Hi Guys,
I have tried contacting sales but it looks like they don’t answer questions about Student or Home edition!
I am looking to get MatLab as a home user in order to learn MatLab/Simulink.
I am not doing a course at an Educational establishment but rather studying at home in my own time.
I’m interested in the following:
MatLab, Simulink, Audio System Toolbox, Signal Processing Toolbox, DSP System Toolbox, MatLab Coder and Simulink Coder
The problem is I cannot seem to add MatLab Coder and Simulink Coder to the order, they are missing in the Home edition although they are there in the student edition.
Is it not possible for me to get the Coders?
Cheers
AndyHi Guys,
I have tried contacting sales but it looks like they don’t answer questions about Student or Home edition!
I am looking to get MatLab as a home user in order to learn MatLab/Simulink.
I am not doing a course at an Educational establishment but rather studying at home in my own time.
I’m interested in the following:
MatLab, Simulink, Audio System Toolbox, Signal Processing Toolbox, DSP System Toolbox, MatLab Coder and Simulink Coder
The problem is I cannot seem to add MatLab Coder and Simulink Coder to the order, they are missing in the Home edition although they are there in the student edition.
Is it not possible for me to get the Coders?
Cheers
Andy Hi Guys,
I have tried contacting sales but it looks like they don’t answer questions about Student or Home edition!
I am looking to get MatLab as a home user in order to learn MatLab/Simulink.
I am not doing a course at an Educational establishment but rather studying at home in my own time.
I’m interested in the following:
MatLab, Simulink, Audio System Toolbox, Signal Processing Toolbox, DSP System Toolbox, MatLab Coder and Simulink Coder
The problem is I cannot seem to add MatLab Coder and Simulink Coder to the order, they are missing in the Home edition although they are there in the student edition.
Is it not possible for me to get the Coders?
Cheers
Andy matlab home, home license MATLAB Answers — New Questions
matlab RIS MBM Communication systems
Matlab完成RIS-MBM的相关系统设计,使系统在发射天线数N=20,接收天线数Nr=3,,RF镜个数为2的条件下,信噪比SNR=15dB时,使用QPSK调制条件下BER低于10-2。Matlab完成RIS-MBM的相关系统设计,使系统在发射天线数N=20,接收天线数Nr=3,,RF镜个数为2的条件下,信噪比SNR=15dB时,使用QPSK调制条件下BER低于10-2。 Matlab完成RIS-MBM的相关系统设计,使系统在发射天线数N=20,接收天线数Nr=3,,RF镜个数为2的条件下,信噪比SNR=15dB时,使用QPSK调制条件下BER低于10-2。 matlab MATLAB Answers — New Questions
pso algoritm in pv
i want to write pso algorithm for estimate . please help me about cost functioni want to write pso algorithm for estimate . please help me about cost function i want to write pso algorithm for estimate . please help me about cost function pso, pv MATLAB Answers — New Questions
passing variable from matlab to python using pyrunfile
Hi all,
I’ve been looking to find an answer for the following question;
I have a variable in my workspace that I need to use in a python script, which for some reason is more challenging in than expected. I have been looking around the Mathlab documentation which unfortunately has not been helpfull for me.
Any help to solve this would be greatly appreciated!
Sample code which I simplified for this message.
% Matlab: This is the line that calls function.py using var01, var02 and var03.
apply_cnn_command = pyrunfile ("function.py ‘" + var01 + var02 + var03+"’")
% In Python I am using the current code:
import sys
import os
import matlab.engine
eng= matlab.engine.start_matlab()
varpy01 = eng.workspace[var01]
% The error message: "Python Error: NameError: name var01 is not defined.Hi all,
I’ve been looking to find an answer for the following question;
I have a variable in my workspace that I need to use in a python script, which for some reason is more challenging in than expected. I have been looking around the Mathlab documentation which unfortunately has not been helpfull for me.
Any help to solve this would be greatly appreciated!
Sample code which I simplified for this message.
% Matlab: This is the line that calls function.py using var01, var02 and var03.
apply_cnn_command = pyrunfile ("function.py ‘" + var01 + var02 + var03+"’")
% In Python I am using the current code:
import sys
import os
import matlab.engine
eng= matlab.engine.start_matlab()
varpy01 = eng.workspace[var01]
% The error message: "Python Error: NameError: name var01 is not defined. Hi all,
I’ve been looking to find an answer for the following question;
I have a variable in my workspace that I need to use in a python script, which for some reason is more challenging in than expected. I have been looking around the Mathlab documentation which unfortunately has not been helpfull for me.
Any help to solve this would be greatly appreciated!
Sample code which I simplified for this message.
% Matlab: This is the line that calls function.py using var01, var02 and var03.
apply_cnn_command = pyrunfile ("function.py ‘" + var01 + var02 + var03+"’")
% In Python I am using the current code:
import sys
import os
import matlab.engine
eng= matlab.engine.start_matlab()
varpy01 = eng.workspace[var01]
% The error message: "Python Error: NameError: name var01 is not defined. pyrunfile, syntax, workspace variables MATLAB Answers — New Questions
Problem with converted TIFF image
I converted a png to tiff unsigned int16. The image resulting is completely black.
This a apart of my code:
% Convert to unsigned 16-bit integer
img_uint16 = uint16(double(img) / double(intmax(‘uint16’)) * 65535);
% Save the image as TIFF
tiffPath = fullfile(subfolderPath, [currentImage(1:end-4), ‘_converted.tif’]);
imwrite(img_uint16, tiffPath);I converted a png to tiff unsigned int16. The image resulting is completely black.
This a apart of my code:
% Convert to unsigned 16-bit integer
img_uint16 = uint16(double(img) / double(intmax(‘uint16’)) * 65535);
% Save the image as TIFF
tiffPath = fullfile(subfolderPath, [currentImage(1:end-4), ‘_converted.tif’]);
imwrite(img_uint16, tiffPath); I converted a png to tiff unsigned int16. The image resulting is completely black.
This a apart of my code:
% Convert to unsigned 16-bit integer
img_uint16 = uint16(double(img) / double(intmax(‘uint16’)) * 65535);
% Save the image as TIFF
tiffPath = fullfile(subfolderPath, [currentImage(1:end-4), ‘_converted.tif’]);
imwrite(img_uint16, tiffPath); tiff MATLAB Answers — New Questions
LS method in idetification ,control engeeniring
Hi,how are you?I have problem with that code,is about LS method with ARX model in identification,I should use xcorr command for EPS variable to compare it that how much is near to white noise ?Hi,how are you?I have problem with that code,is about LS method with ARX model in identification,I should use xcorr command for EPS variable to compare it that how much is near to white noise ? Hi,how are you?I have problem with that code,is about LS method with ARX model in identification,I should use xcorr command for EPS variable to compare it that how much is near to white noise ? control engeeniring MATLAB Answers — New Questions
FEM programs to solve ODE
I can’t understand the ‘feasmbl2’. Is that ok to map the matrix and vectors of different sizes?
It happens the error context as the following below:
error: ff(_,2): out of bound 1 (dimensions are 6×1)
error: called from
feasmbl2 at line 30 column 18
Ex351 at line 83 column 8
Here is my whole codes from the textbook, but it can’t run well.
%—————————————————————————-
% EX3.5.1
% to solve the ordinary differential equation given as
% a u” + b u’ + c u = 1, 0 < x < 1
% u(0) = 0 and u(1) = 0
% using 5 linear elements
%
% Variable descriptions
% k = element matrix
% f = element vector
% kk = system matrix
% ff = system vector
% index = a vector containing system dofs associated with each element
% bcdof = a vector containing dofs associated with boundary conditions
% bcval = a vector containing boundary condition values associated with
% the dofs in ‘bcdof’
%—————————————————————————-
%————————————
% input data for control parameters
%————————————
clear
nel=5; % number of elements
nnel=2; % number of nodes per element
ndof=1; % number of dofs per node
nnode=6; % total number of nodes in system
sdof=nnode*ndof; % total system dofs
%—————————————–
% input data for nodal coordinate values
%—————————————–
gcoord(1)=0.0; gcoord(2)=0.2; gcoord(3)=0.4; gcoord(4)=0.6;
gcoord(5)=0.8; gcoord(6)=1.0;
%—————————————————–
% input data for nodal connectivity for each element
%—————————————————–
nodes(1,1)=1; nodes(1,2)=2; nodes(2,1)=2; nodes(2,2)=3;
nodes(3,1)=3; nodes(3,2)=4; nodes(4,1)=4; nodes(4,2)=5;
nodes(5,1)=5; nodes(5,2)=6;
%—————————————–
% input data for coefficients of the ODE
%—————————————–
acoef=1; % coefficient ‘a’ of the diff eqn
bcoef=-3; % coefficient ‘b’ of the diff eqn
ccoef=2; % coefficient ‘c’ of the diff eqn
%————————————-
% input data for boundary conditions
%————————————-
bcdof(1)=1; % first node is constrained
bcval(1)=0; % whose described value is 0
bcdof(2)=6; % 6th node is constrained
bcval(2)=0; % whose described value is 0
%—————————————–
% initialization of matrices and vectors
%—————————————–
ff=zeros(sdof,1); % initialization of system force vector
kk=zeros(sdof,sdof); % initialization of system matrix
index=zeros(nnel*ndof,1); % initialization of index vector
%—————————————————————–
% computation of element matrices and vectors and their assembly
%—————————————————————–
for iel=1:nel % loop for the total number of elements
nl=nodes(iel,1); nr=nodes(iel,2); % extract nodes for (iel)-th element
xl=gcoord(nl); xr=gcoord(nr);% extract nodal coord values for the element
eleng=xr-xl; % element length
index=feeldof1(iel,nnel,ndof);% extract system dofs associated with element
k=feode2l(acoef,bcoef,ccoef,eleng); % compute element matrix
f=fef1l(xl,xr); % compute element vector
[kk,ff]=feasmbl2(kk,ff,k,f,index); % assemble element matrices and vectors
end
%—————————–
% apply boundary conditions
%————————-
[kk,ff]=feaplyc2(kk,ff,bcdof,bcval);
%—————————-
% solve the matrix equation
%—————————-
fsol=kkff;
%———————
% analytical solution
%———————
c1=0.5/exp(1);
c2=-0.5*(1+1/exp(1));
for i=1:nnode
x=gcoord(i);
esol(i)=c1*exp(2*x)+c2*exp(x)+1/2;
end
%————————————
% print both exact and fem solutions
%————————————
num=1:1:sdof;
store=[num’ fsol esol’]
%—————————————————————
Thanks a lot for helping to debug the code!!I can’t understand the ‘feasmbl2’. Is that ok to map the matrix and vectors of different sizes?
It happens the error context as the following below:
error: ff(_,2): out of bound 1 (dimensions are 6×1)
error: called from
feasmbl2 at line 30 column 18
Ex351 at line 83 column 8
Here is my whole codes from the textbook, but it can’t run well.
%—————————————————————————-
% EX3.5.1
% to solve the ordinary differential equation given as
% a u” + b u’ + c u = 1, 0 < x < 1
% u(0) = 0 and u(1) = 0
% using 5 linear elements
%
% Variable descriptions
% k = element matrix
% f = element vector
% kk = system matrix
% ff = system vector
% index = a vector containing system dofs associated with each element
% bcdof = a vector containing dofs associated with boundary conditions
% bcval = a vector containing boundary condition values associated with
% the dofs in ‘bcdof’
%—————————————————————————-
%————————————
% input data for control parameters
%————————————
clear
nel=5; % number of elements
nnel=2; % number of nodes per element
ndof=1; % number of dofs per node
nnode=6; % total number of nodes in system
sdof=nnode*ndof; % total system dofs
%—————————————–
% input data for nodal coordinate values
%—————————————–
gcoord(1)=0.0; gcoord(2)=0.2; gcoord(3)=0.4; gcoord(4)=0.6;
gcoord(5)=0.8; gcoord(6)=1.0;
%—————————————————–
% input data for nodal connectivity for each element
%—————————————————–
nodes(1,1)=1; nodes(1,2)=2; nodes(2,1)=2; nodes(2,2)=3;
nodes(3,1)=3; nodes(3,2)=4; nodes(4,1)=4; nodes(4,2)=5;
nodes(5,1)=5; nodes(5,2)=6;
%—————————————–
% input data for coefficients of the ODE
%—————————————–
acoef=1; % coefficient ‘a’ of the diff eqn
bcoef=-3; % coefficient ‘b’ of the diff eqn
ccoef=2; % coefficient ‘c’ of the diff eqn
%————————————-
% input data for boundary conditions
%————————————-
bcdof(1)=1; % first node is constrained
bcval(1)=0; % whose described value is 0
bcdof(2)=6; % 6th node is constrained
bcval(2)=0; % whose described value is 0
%—————————————–
% initialization of matrices and vectors
%—————————————–
ff=zeros(sdof,1); % initialization of system force vector
kk=zeros(sdof,sdof); % initialization of system matrix
index=zeros(nnel*ndof,1); % initialization of index vector
%—————————————————————–
% computation of element matrices and vectors and their assembly
%—————————————————————–
for iel=1:nel % loop for the total number of elements
nl=nodes(iel,1); nr=nodes(iel,2); % extract nodes for (iel)-th element
xl=gcoord(nl); xr=gcoord(nr);% extract nodal coord values for the element
eleng=xr-xl; % element length
index=feeldof1(iel,nnel,ndof);% extract system dofs associated with element
k=feode2l(acoef,bcoef,ccoef,eleng); % compute element matrix
f=fef1l(xl,xr); % compute element vector
[kk,ff]=feasmbl2(kk,ff,k,f,index); % assemble element matrices and vectors
end
%—————————–
% apply boundary conditions
%————————-
[kk,ff]=feaplyc2(kk,ff,bcdof,bcval);
%—————————-
% solve the matrix equation
%—————————-
fsol=kkff;
%———————
% analytical solution
%———————
c1=0.5/exp(1);
c2=-0.5*(1+1/exp(1));
for i=1:nnode
x=gcoord(i);
esol(i)=c1*exp(2*x)+c2*exp(x)+1/2;
end
%————————————
% print both exact and fem solutions
%————————————
num=1:1:sdof;
store=[num’ fsol esol’]
%—————————————————————
Thanks a lot for helping to debug the code!! I can’t understand the ‘feasmbl2’. Is that ok to map the matrix and vectors of different sizes?
It happens the error context as the following below:
error: ff(_,2): out of bound 1 (dimensions are 6×1)
error: called from
feasmbl2 at line 30 column 18
Ex351 at line 83 column 8
Here is my whole codes from the textbook, but it can’t run well.
%—————————————————————————-
% EX3.5.1
% to solve the ordinary differential equation given as
% a u” + b u’ + c u = 1, 0 < x < 1
% u(0) = 0 and u(1) = 0
% using 5 linear elements
%
% Variable descriptions
% k = element matrix
% f = element vector
% kk = system matrix
% ff = system vector
% index = a vector containing system dofs associated with each element
% bcdof = a vector containing dofs associated with boundary conditions
% bcval = a vector containing boundary condition values associated with
% the dofs in ‘bcdof’
%—————————————————————————-
%————————————
% input data for control parameters
%————————————
clear
nel=5; % number of elements
nnel=2; % number of nodes per element
ndof=1; % number of dofs per node
nnode=6; % total number of nodes in system
sdof=nnode*ndof; % total system dofs
%—————————————–
% input data for nodal coordinate values
%—————————————–
gcoord(1)=0.0; gcoord(2)=0.2; gcoord(3)=0.4; gcoord(4)=0.6;
gcoord(5)=0.8; gcoord(6)=1.0;
%—————————————————–
% input data for nodal connectivity for each element
%—————————————————–
nodes(1,1)=1; nodes(1,2)=2; nodes(2,1)=2; nodes(2,2)=3;
nodes(3,1)=3; nodes(3,2)=4; nodes(4,1)=4; nodes(4,2)=5;
nodes(5,1)=5; nodes(5,2)=6;
%—————————————–
% input data for coefficients of the ODE
%—————————————–
acoef=1; % coefficient ‘a’ of the diff eqn
bcoef=-3; % coefficient ‘b’ of the diff eqn
ccoef=2; % coefficient ‘c’ of the diff eqn
%————————————-
% input data for boundary conditions
%————————————-
bcdof(1)=1; % first node is constrained
bcval(1)=0; % whose described value is 0
bcdof(2)=6; % 6th node is constrained
bcval(2)=0; % whose described value is 0
%—————————————–
% initialization of matrices and vectors
%—————————————–
ff=zeros(sdof,1); % initialization of system force vector
kk=zeros(sdof,sdof); % initialization of system matrix
index=zeros(nnel*ndof,1); % initialization of index vector
%—————————————————————–
% computation of element matrices and vectors and their assembly
%—————————————————————–
for iel=1:nel % loop for the total number of elements
nl=nodes(iel,1); nr=nodes(iel,2); % extract nodes for (iel)-th element
xl=gcoord(nl); xr=gcoord(nr);% extract nodal coord values for the element
eleng=xr-xl; % element length
index=feeldof1(iel,nnel,ndof);% extract system dofs associated with element
k=feode2l(acoef,bcoef,ccoef,eleng); % compute element matrix
f=fef1l(xl,xr); % compute element vector
[kk,ff]=feasmbl2(kk,ff,k,f,index); % assemble element matrices and vectors
end
%—————————–
% apply boundary conditions
%————————-
[kk,ff]=feaplyc2(kk,ff,bcdof,bcval);
%—————————-
% solve the matrix equation
%—————————-
fsol=kkff;
%———————
% analytical solution
%———————
c1=0.5/exp(1);
c2=-0.5*(1+1/exp(1));
for i=1:nnode
x=gcoord(i);
esol(i)=c1*exp(2*x)+c2*exp(x)+1/2;
end
%————————————
% print both exact and fem solutions
%————————————
num=1:1:sdof;
store=[num’ fsol esol’]
%—————————————————————
Thanks a lot for helping to debug the code!! fem, feasmbl2, out of bound, ode MATLAB Answers — New Questions