Writing and Reading Simple Text Data (MATLAB)

Working with text data is convenient as we can share our data with other systems without giving them the trouble to install a particular software to read our data.

Let’s see how we can read and write a text data using MATLAB.

%% Reading the data in text file

fileID = fopen('nums.txt','r'); %opeing the file in reading mode

formatSpec = '%d %f'; %defining the format of the data

sizeA = [2 Inf]; %defining the size of the data

A = fscanf(fileID,formatSpec,sizeA); %reading the data using fscanf function

fclose(fileID); %closing the file

A' %displaying the data

Writing Text Data

clear; close all; clc;

%% Defining the variables

x = 1:10;

y = [x;10*rand(1,10)]; %randomly generated 2nd column

z=['x-values' '  ' 'random-values']; %writing the header

%% Writing the data in text file

fileID = fopen('nums.txt','w'); %opens a text file in writing mode

fprintf(fileID,'%s %s %s\n',z);

fprintf(fileID,'\n');

fprintf(fileID,' %d         %4.4f\n',y); %write the data y in the file.

%The first column is written in digit format and second column in floating point format.

fclose(fileID); %close the file

MATLAB Scripts:

Writing text data

Reading Data from Microsoft Excel (MATLAB)

The MATLAB is very easy to use language. This is one of the main reason why it is getting popular day by day.

Here, we show how we can read data from the Microsoft excel into the MATLAB. Once we can read the data into the MATLAB, we can do other computations and operations quite easily.

Sample data file (global_temp.xlsx)

Saving data in Mat or text format for easier and faster use in later programs

Constrained Least Square Fitting

Let us try to understand the least square fitting method using an example:

Example

Seismic Resolution

Seismic resolution and fidelity are the two important measures of the quality of the seismic record and the seismic images. Seismic resolution quantifies the level of precision, such as the finest size of the subsurface objects detectable by the seismic data whereas the seismic fidelity quantifies the truthfulness such as the genuineness of the data or the level to which the imaged target position matches its true subsurface position.

Let us try to understand this by making a synthetic data and doing the analysis over it.

Seismic Resolution Analysis: 1

Estimation of the Width of the Fresnel’s Zone

The Fresnel’s resolution quantifies the resolvability of seismic wave perpendicular to the direction of wave propagation. Fresnel’s resolution is defined as the width of the first Fresnel’s zone due to interference of the spherical waves from the source and from the receiver.

Seismic Resolution Analysis: 2

Fidelity

To investigate the fidelity of our data, let us consider the technique of resampling. For our case we consider the method of “bootstrapping”. Bootstrapping basically relies on random sampling with replacement. The other popular method for resampling is “jackknifing” which predates “bootstrapping”. The jackknife estimator of a parameter is found by systematically leaving out each observation from a dataset and calculating the estimate and then finding the average of these calculations.

The principle behind “bootstrapping” is that a dataset is taken, the total dataset is divided into two by randomly sampling with replacement. The newly sampled data are now used to invert for the model using some kernel function. If the two models correlate high enough then we can say that the prominent features in the model come from consistent signals in the data.

We don’t have the data set to make the velocity model. So instead, we can take random Gaussian distribution data and play with it.

Let’s pose the null hypothesis that the two sets of data come from the same probability distribution (not necessarily Gaussian). Under the null hypothesis, the two sets of data are interchangeable, so if we aggregate the data points and randomly divide the data points into two sets, then the results should be comparable to the results obtained with the original data. So, the strategy is to generate random datasets, with replacement (bootstrapping), compute difference in means (or difference in medians or any other reliable statistic), and then compare the resulting values to the statistic computed from the original data.

Seismic Resolution Analysis: 3

Analytical signal and Hilbert Transform

For a real time series x(t), its analytic signal x(t) is defined as

x(t) = x(t) – iH[x(t)]

Let us consider an example of a monochromatic signal 𝑥(𝑡) = 5 sin(10𝑡 + 3).

Figure 1

Now, let us consider a more complex function

x(t) = 1*sin(2𝜋10𝑡 + 0.3) + 2*sin(2𝜋20𝑡 + 0.2) + 3*sin(2𝜋30𝑡 + 0.4).

Figure 2

We can clearly observe that the Hilbert transform estimates the instantaneous frequency of a signal for monocomponent signals only.

Matlab Codes:

clear; close all; clc

fs = 1e4;

t = 0:1/fs:1;

f1=10;

f2=20;

f3=30;

x=1*sin(2*pi*f1*t+0.3)+2*sin(2*pi*f2*t+0.2)+3*sin(2*pi*f3*t+0.4);

%x=5*sin(2*pi*10*t+3);

y = hilbert(x);

figure(1)

plot(t,real(y),t,imag(y))

%xlim([0.01 0.03])

legend('real','imaginary')

title('Hilbert Function')

figure(2)

subplot(3,1,1)

env=abs(y);

plot(t,x)

xlabel('Time')

title('Envelope')

hold on

plot(t,env)

legend('original','envelope')

subplot(312)

instph=fs/(2*pi)*unwrap(angle(y));

plot(t,instph)

xlabel('Time')

grid on

title('Instantaneous Phase')

subplot(313)

instfreq = fs/(2*pi)*diff(unwrap(angle(y)));

plot(t(2:end),instfreq)

xlabel('Time')

ylabel('Hz')

grid on

title('Instantaneous Frequency')

Application of Hilbert Transform:

Figure 3

Aliasing and Nyquist Condition

Let us consider a simple sinusoidal signal

x(t) = 5sin(10t + 3)

This is a signal with an amplitude of 5, the angular frequency of 10 and phase of 3. The Nyquist condition demands that the angular frequency, mod(𝜔) ≤ pi/delta t.

For the above figure 1, we took 𝜔 to be 10 and pi/delta t is 99.5, so the Nyquist condition satisfies and we do not see the aliasing. Here the sampling rate, Δ𝑡 is 0.0316.

Let us see what happens when we take the pi/delta t to be less than 𝜔.

In this case, we took the sampling rate, Δ𝑡 to be 0.3307. Here the pi/delta t is 9.5000 which is less than 𝜔. Hence the Nyquist condition is not satisfied and there is a clear aliasing. The sinusoidal signal does not appear to be sinusoidal.

Figures