How can I get the Stockwell transform of an earthquake data in MATLAB?

Dear, I was wondering does MATLAB have any tool to calculate the Stockwell transform ? If it exists in MATLAB library, can you please provide me a script code so I can get S-tranform of an earthquake data. Thanks for your help. Regards

Answer 1

[dpb]

See if one of the following File Exchange submittals suits...

<Stockwell Transform (S-Transform)>

<FFT-fast S-transforms: DOST, DCST, DOST2 and DCST2>

<Time frequency distribution of a signal using S-transform (stockwell transform)>

[Adnan]

Thanks for that, yet each one has some bugs. I got some errors from their codes. Can we not find a correct code ? Can you please provide me a correct one ? Thanks again for your time sir.

[dpb]

The code for the last above does run, albeit it is pretty lacking in style and formatting that could make it much easier to read. I did take the time to insert the missing end statements at the end of the various functions; there should be one for st at about line 145 (just after a return statement) that isn't needed.

Then another for function strans about 235; there's a comment about "end strans function".

About line 263 for function g_window

About line 305 for function check_input (comment about "bloody odd how you don't end function" shows lack of MATLAB expertise by the original author)

Last for function st_help at the very end of the file.

Not mandatory, but a lot cleaner code that way...the one thing is if add one, must add them all; can't mix metaphors in style within on m file.

Answer 2

[dpb]

What errors?

I have none of these downloaded and don't know of any specific others...and, personally, am not familiar with the Stockwell transform.

Well, I downloaded the last of the above but used the original st function directly rather than the poster's version that uses imcomplement on the output as I didn't see the purpose.

I then used the quake.mat file distributed with MATLAB and ran the E-W acceleration through it ---

>> EQ=load('quake');
>> [ST,t,f]=st(EQ.e,0,100,1/200);
 
Minfreq = 0
Maxfreq = 100
Sampling Rate (time   domain) = 5.000000e-03
Sampling Rate (freq.  domain) = 1
The length of the timeseries is 10001 points
 
The number of frequency voices is 101
Calculating analytic signal (using Hilbert transform)
Estimated time is 0.078174Calculating S transform...
Finished Calculation
>> imagesc(t,f,abs(ST))
>> xlabel('Frequency, Hz')
>> ylabel('Time (s)')
>> title('Oct17, 1989 Loma Prieta earthquake E-W ground motion')
>> 

with the result attached...looks as though it worked to me, albeit I've never used the transform so don't really know anything about it other than it is another of the time-frequency separation tools.

ADDENDUM

Above with the suggested change to use

analytic_signal=isreal(timeseries);

instead of the hardcoded FALSE and the correct Nyquist frequency for the actual input timeseries. You'll observe the image is the same other than the y axis labels...

[Adnan]

Thanks for that, yet I could not run again. I was just wondering, if you do not mind, can you please share that code with me? Thanks again. It can be in zip file or normal scrip in this answer box.

[dpb]

It's the code at the FEX with only cosmetic fixes; post the error message and the code you ran that gave an error....

I do see one thing I would change; at the beginning Mr. Stockwell set the analytic_signal flag variable to FALSE; I would change that to use

analytic_signal=isreal(timeseries);

so it will be based on the actual input data array.

There is one typo in the input for the specific Loma Prieta earthquake data I used for example -- that is the sampling rate was 200 Hz (time 1/200) so the max frequency should have been 100, not 500. That only changes the axis scaling on the output plot; returned frequency vector.

[dpb]

See <Stockwell Transform by Robert Stockwell MATLAB demo>

[Adnan]

Thanks, @dpb for that, yet it does not work. There is only one ".m" file in your file, and it is a function. Can you please add the other file if you do not mind? I do appreciate your help.

[dpb]

It DOES work and there is only the one file. Call it as the help and my code example shows.

Of course, the function is now named stockwell_transform_st to distinguish this FEX submittal from previous ones.

EQ=load('quake');
[ST,t,f]=stockwell_transform_st(EQ.e,0,100,1/200);

[Adnan]

Thanks again for your valuble time dr. @dpb, I have earthquake data; therefore, I run this script (below; Script No.1) apart from your code. I got this error (given below). Can you please fix that for me?

% error obtained from MATLAB
error: "Dot indexing is not supported for variables of this type"

% runing part (Script No.1) 
load DATA.txt
EQ=DATA(:,1);
[ST,t,f]=stockwell_transform_st(EQ.e,0,100,1/200);

[dpb]

You have to pay attention to what you're doing and use your own data correctly...

The EQ struct is the result of my having loaded the MATLAB .mat file -- it contains data arrays inside it that are accelerations. Look at the documentation example in MATLAB to see that -- it's an example for using a timetable but that's immaterial here; as it will show, there are three variables e,n,v in the .mat file, one for the E-W, N-S and Vertical motion accelerometers. Thus "EQ.e" above is returning the EW accelerometer data.

You have, apparently, some data in a text file; only you know what that is, but you can't reference an array with a reference to a struct.

And, you need to make the other input parameters match whatever the sampling characteristics were for the data you do have; the MATLAB example data were sampled at 200 Hz; there's no guarantee your data were also, but we can't/don't know that.

% runing part (Script No.1) 
load DATA.txt
EQ=DATA(:,1);
minfreq=???             % what do you want for the minimum analysis frequency
maxfreq=???             % and the max (Nyquist is as high as you can go)??
samplingrate=???        % what was the sampling rate of the signal (delta-t)

Once you have that defined, THEN you can call

[ST,t,f]=stockwell_transform_st(EQ,minfreq,maxfreq,samplingrate);

You have to do some work on your end, here, too...things don't just happen by magic. You need to read the help for the function and probably study how it is put together at least some so you have some idea of what you're using, not just as a "black box"

[Adnan]

Thank you very much for that sir @dpb. Unfortunately, the results for the earthquake record do not make sense (I compared them with the published paper results); perhaps this code is suitable for signal-related matters.

I was just wondering, if you don't mind, I have attached the Stockwell equations here (a pdf file). Could you please assist me in generating this equation in Matlab for the earthquake record (DATA.txt: the first column represents time and the second column represents acceleration)? I extracted the equations (equations 1 to 8) from a paper, and I am confident this equation will assist other researchers in accessing similar information in the MATLAB community file, and it will likely be downloaded numerous times.

I look forward to receiving your reply, sir.

Regards,

P.S. the paper name that I extracted the equation is "On the use of Stockwell transform in structural dynamic analysis"

[dpb]

The code you have was written by Stockwell himself...don't believe you'll get any closer to the oracle than that, and having looked at that code, it is an implementation of a segmented FFT over the time trace.

As the help/documentation with the file notes, there are user-adjustable parameters buried inside the function that can be adjusted; you may need to compare what the author of the particular paper chose for such details and adjust those and, potentially, see if there are differences in normalization choices and/or whether the acceleration data are in particular units (g, m/s^2, V, ...?).

Attach the file and supporting information including the code you used to run it; we can't see a thing on your desktop here.

"... DATA.txt: the first column represents time and the second column represents acceleration"

In your prior posting

load DATA.txt
EQ=DATA(:,1);
[ST,t,f]=stockwell_transform_st(EQ.e,0,100,1/200);

you passed the first column, not the second...did you fix that? Again, w/o seeing what you actually did and knowing the particulars of the dataset, we're blind.

Finally, regarding "...Could you please assist me in generating this equation in Matlab..?"

The Answers forum is not a free personal consulting platform, my standard consulting rates are $100/hr if you want more.

[Adnan]

I see, Thanks for your reply sir. It is a lot to me. Here it is my code and the data that I used. I added the earthqauke record (ImperialValley.txt- the record is attached below at the bottom of this email). Once I run this code. The output should be like around Image_No_1.PNG (attached here). Thanks again sir for your time.

Output part:

%% Main RUN 
load ImperialValley.txt; %% the data is attached below (at the bottom of this email) 
EQ=Imperial(:,2);
time=Imperial(:,1);
dt = mean(diff(time));
minfreq=-10;           
maxfreq=10    ;        
samplingrate=1/dt ;    

[ST,t,f]=stockwell_transform_st(EQ,minfreq,maxfreq,samplingrate);
figure;
imagesc(f,t,abs(ST))
xlabel('Frequency, Hz')
ylabel('Time (s)')

function part:

%% stockwell_transform_st

Removed FEX code; not comfortable posting in forum -- dpb

[dpb]

%% Main RUN 
load ImperialValley.txt; %% the data is attached below (at the bottom of this email) 
EQ=Imperial(:,2);
time=Imperial(:,1);
dt = mean(diff(time));
minfreq=-10;           
maxfreq=10    ;        
samplingrate=1/dt ;    

[ST,t,f]=stockwell_transform_st(EQ,minfreq,maxfreq,samplingrate);

I'd probably do the initial something like

tIV=readtable('ImperialValley.txt'); % tables are convenient for viewing as well as using...
tIV.Properties.VariableNames={'Time','Accel'};  % add convenient names
head(tIV)                                       % take a look at what we have...
dt=mean(diff(tIV.Time));                        % now use the table variables
plot(tIV.Time,tIV.Accel)                        % look at the time trace
...

That's just convenient way to have the data; has no bearing on the actual calculations other than what choice make in how to read/handle it...

As to the transform, as noted before I have never used the Stockman transform (hadn't actually even heard of it before), so I'm certainly not well versed in it or its application and all I can do is presume Stockman knew how to implement what he intended.

That said, there are some things in the code that seem peculiar but I don't have the time to try to research the topic sufficiently to fully understand -- the foremost of those is his use of the FREQSAMPLINGRATE input and resultant calculation of the returned frequency vector; you'll need to dig into that and understand what is being done there.

minfreq=-10;

There is no such thing as a negative frequency here; there is in a two-sided FFT, but not here. It appears you've taken the exponent of the lower limit of the log range of the sample figure here instead; if you note the output of the function, it caught the negative Minfreq and set it to zero.

maxfreq=10    ;        
samplingrate=1/dt ;    

The sampling rate input is documented in the usage as

"SAMPLINGRATE - the time interval between successive data points"

so it should be dt, not 1/dt. Changing that will make the time axis of the plot look as you expect. It looks like you again copied my prior code where I used 1/Fs, but Fs was the sampling RATE of the signal, 200 Hz, so 1/FS --> dt.

There is an issue in the way Mr/Dr? Stockwell wrote this code in that one cannot use fractional frequencies such as 0.1 for minFreq or maxFreq; the code uses an internal colon operation between the two that breaks if the difference isn't an integer.

I suspect one could figure out and modify this code base to handle the binning more elegantly and with investigation of the internal width factor on the gaussian window to get what appears to be finer resolution in the time domain in the reference; as noted above, those details between this code and that author's are all pieces that will take researching to be able to duplicate one from the other.

As noted, there are other implementations of what appear to be the same idea on the File Exchange, there's also a C version of the algorithm I found in a google search that can be turned into a mex file and called from MATLAB.

I'm sorry, but I simply don't have the time to spare to investigate more fully; if you need something like this and don't feel you have the programming skills to do so yourself, you might try contacting the author of the other paper and see if they have supporting code available.

ERRATUM:

The output display of the version I uploaded will be much more readable if you will search for each line containing fprintf and add a \n to the end of the format string; one cleanup I did was to use the MATLAB suggestion to fix the use of the disp(sprintf(...)) idiom automatically; the problem is their builtin fix doesn't add the trailing newline into the format string that disp() inserts automagically so those lines run on instead of showing as separate lines...at some point I'll upload a new version with that correction, but there are <10 lines and it's a trivial edit to make.

Good luck...

[Adnan]

thanks for that sir @dpb. I will consider your recommendation!.