Matlab oppgave 1

Knut Sørsdal GEO 4280 Matlab Exc. 1.

for n = 1:16

x(n)=randn

end

stem(x)

%Compute DFT with fft

Rex=real(fft(x))

Imx=imag(fft(x))

absx=abs(fft(x))

Phasex=angle(fft(x))

y=0:15

subplot(2,2,1), stem(y,Rex)

title Re(DFT(x(n)))

xlabel frequency.points

subplot(2,2,2), stem(y,Imx)

title Im(DFT(x(n)))

xlabel frequency.points

subplot(2,2,3), stem(y,absx)

title Magnitude(DFT(x(n)))

xlabel frequency.points

subplot(2,2,4), stem(y,Phasex)

title Phase(DFT(x(n)))

xlabel frequency.points

C-component is the value for zero-frequency. In this solution I have a serie of L=16 points. Frequency point L/2 =8 is called the Nyquist-frequency. The DFT have 14 regular points + DC + Nyquist frequency with symmetry around Nyquistfrequency for Re(DFT(x(n)) and antisymmetry for Im(DFT(x(n). To see the symmetry (and antisymmetry) we would have to introduce negative frequencies (that is not done on the plot.). Then we would have seen values above L/2 repeated from –L/2 up to DC.

Discussion of Nyquist-frequency: We need at least two sample points for a wave-cycle to avoid aliasing for that frequency, and that is called Nyquist-criterion. 16 points divided 2 gives Nyquist-frequency fn=8 points. If we had defined the length of x(n) as time and x(n) as a time-series we would have a Nyquist frequency 8 Hz. The important point is that the specter we will see around DC (included negative frequencies), will repeat itself from Nyquist-frequency and above.

c) The same calculations for N=17. First plot is x(n)

We can use same DFT calculations for x(n) of length N=17 as in b) and get the result:

We can have the same discussion as on b) and get very similar conclusions and we have a DC-component. But L/2 gives not a Nyquist-frequency when L=17 (at least we can not see it). We can conclude that we cannot see the Nyquist-frequency when the length of the series x(n) is an odd number.
Exc. 2

clear

l=21

box=boxcar(l)

black=blackman(l)

x=1:21

plot(x,box,x,black)

title rectangular(blue)/blackman(green)

a=4

m=(l-1)/2

%Cauchy-window

for n=1:21

w(n)=(m*m)/(m*m + a*a*(n-m)*(n-m))

end

plot(x,box,x,black,x,w)

title time.windows

legend(‘rectangular’,’blackman’,’Cauchy’)

%Compute DFT with fft padding to 101

absblack=abs(fft(black,101))

absbox=abs(fft(box,101))

absw=abs(fft(w,101))

%plot

x=1:25

subplot(2,1,1), plot(x,absbox(1:25),x,absblack(1:25),x,absw(1:25))

title Magnitude

xlabel frequency

legend('rectangular','blackman','Cauchy')

x=1:21

subplot(2,1,2), plot(x,box,x,black,x,w)

title time-windows

xlabel time

legend('rectangular','blackman','Cauchy')

%Cauchy for different values of a

for n=1:21

w1(n)=( m*m)/(m*m + 1*1*(n-m)*(n-m))

w2(n)=( m*m)/(m*m + 2*2*(n-m)*(n-m))

w3(n)=( m*m)/(m*m + 3*3*(n-m)*(n-m))

w4(n)=( m*m)/(m*m + 4*4*(n-m)*(n-m))

w5(n)=( m*m)/(m*m + 5*5*(n-m)*(n-m))

w6(n)=( m*m)/(m*m + 6*6*(n-m)*(n-m))

w7(n)=( m*m)/(m*m + 7*7*(n-m)*(n-m))

w8(n)=( m*m)/(m*m + 8*8*(n-m)*(n-m))

end

x=1:21

plot(x,w1,x,w2, x,w3,x,w4, x,w5,x,w6, x,w7,x,w8)

legend('a=1','a=2','a=3','a=4','a=5','a=6','a=7','a=8')

title 'Cauchy.time-window'

xlabel 'time'

Discussion : Values start on a=1 with a broad time window. When a increases from 1 up to 8 the Cauchy window becomes more narrow.

Esc. 3

Discussion:

Circular convolution gives another solution than linear convolution if the length of the output series is not the same. When they have the same length the solution is also the same. The graph on the right has the same length as the graph on left, and they are similar. The graph in the middle are different and has a different length.
Exc.4 Wiener filter

%This whole program actually runs as is, but all numbers came at the end. They are copied into the program afterwards and colored blue to give an easier reading.

%1)I propose two solutions. The first is with matlabs predefined functions:

s=[6 5 1];

d=[1 0 0];

a = xcorr(s);

a1 = [a(3) a(2) a(1);a(2) a(3) a(2);a(1) a(2) a(3)];

b = xcorr(s,d);

b1 = [b(3);b(2);b(1)];

f = a1\b1;

%filter f with 3 coefficients:

disp(f);

0.1589

-0.1189

0.0518

%The second solution is more inconvenient, but gives same result.

%Autocorrelation of a function s(x) is u_ss=s(z)*s(1/z) that gives a matrix that easily can be found analytically. If m is our prefix, Crosscorrelation u_sd = s₀ for m=0 and 0 elsewhere. After analytic calculations we arrive at these numbers for autocorrelation A and crosscorrelation B:

A = [62 35 6;35 62 35;6 35 62];

%Crosscorrelation s and d gives the solution

B = [6;0;0];

f=A\B;

disp(f);

0.1589

-0.1189

0.0518

%2)

c=conv(s,f);

x=0:4;

subplot(1,5,1), stem(x,c)

disp(c);

%see plot no 1 from left

0.9534

0.0810

-0.1252

0.1398

0.0518

error=((1-c(1))*(1-c(1))) + c(2)*c(2)+c(3)*c(3)+c(4)*c(4)

disp(error)

0.0440

%We have some error energy indicating that we did not get 100 % spiking

%3)

%filter coefficients that we calculated analytic from ex.3 p.11 with γ=1

h=[1/6;-5/36;19/219];

c3=conv(s,h);

error=((1-c3(1))*(1-c3(1))) + c3(2)*c3(2)+c3(3)*c3(3)+c3(4)*c3(4)

subplot(1,5,2), stem(x,c3)

disp(c3)

%see plot no 2 from left

1.0000

-0.0000

-0.0072

0.2949

0.0868

disp(error)

0.0870

% We got more error energy then in 2 and would expect a less correct spiking

%4)

%To get 4 coefficients we must add a zero:

s=[6 5 1 0];

d=[1 0 0 0];

a = xcorr(s);

a1 = [a(4) a(3) a(2) a(1);a(3) a(4) a(3) a(2);a(2) a(3) a(4) a(3); a(1) a(2) a(3) a(4)];

b = xcorr(s,d);

b1 = [b(4);b(3);b(2);b(1)];

f = a1\b1;

%filter f with 4 coefficients:

disp(f);

A1 = [62 35 6 0;35 62 35 6;6 35 62 35;0 6 35 62];

B1 = [6;0;0;0];

f4 = A1\B1;

c4 = conv(s,f4);

error=((1-c4(1))*(1-c4(1))) + c4(2)*c4(2)+c4(3)*c4(3)+c4(4)*c4(4);

disp(f4);

0.1644

-0.1328

0.0761

-0.0301

disp(error);

0.0072

disp(c4);

%see plot no 3 from left

0.9864

0.0250

-0.0434

0.0670

-0.0743

-0.0301

%Now we have much smaller error energy than in 2 and 3, and consequently better spiking.

x=0:6;

subplot(1,5,3), stem(x,c4);

%5)

%We have the same autocorrelation as in 4) with a new input signal:

s= [3 7 2]

% filter coefficients are given from 4):

f4 = [0.1644;-0.1328;0.0761;-0.0301]

c5=conv(s,f4);

x=0:5;

subplot(1,5,4), stem(x,c5)

disp(c5);

%see plot no 4 from left

0.4932

0.7524

-0.3725

0.1768

-0.0585

-0.0602

% we have a bad spiking, much different than from solution under 4.The reason is that our input-signal s is not minimum phase.

%6)

%We have the same autocorrelation as in 4)and 5) and the same s as in 5).We want a two sided filter with 4 coeffisients.

s=[3 7 2 0];

d=[0 0 1 0];

%We still have the same autocorrelation as in 4)(I tested it and it worked)

a1 = [62 35 6 0;35 62 35 6;6 35 62 35;0 6 35 62];

%And we get a new crosscorrelation:

b = xcorr(s,d);

b1 = [b(4);b(3);b(2);b(1)];

f = a1\b1;

disp(f)

-0.0621

0.1761

-0.0520

0.0123

c6 = conv(s,f)

disp(c6)

%see plot no 5 from left

-0.1864

0.0934

0.9527

0.0254

-0.0179

0.0246

subplot(1,5,5), stem(c6)

We now have a considerable better spiking than in 5). The result is better when we delay our filter.

%Matlab lecturer told us to discuss briefly all the graphs depicting the solutions of input series s convolved with filter f.(Total 5)

The best approximation to value 1 for x=0 is on plot number 2 from left, but we could expect less error energy on plot 3 – due to smaller oscillations after spike. That is also the case (from numbers above). We will therefore say that we have best spiking on plot 3. Plot 4 and 5 does not give good spiking because our s is not minimum phase. However plot 5 gives a better spiking than plot 4 because we have delayed the delta-puls in the cross correlation.

Outline for seminary in autumn semester

%A Wienerfilter could be applied on Q-filtered seismic trace to spike minimum phase, zero phase and mix phase Ricker-pulses. I have made these pulses as an appendix to this exc. 4. And hope to apply it correct. I use a matlab packet seislab. I hope to present the data on seminary in autumn 2008. The book of Wang on inverse Q-filtering will also be presented there.

%Here is an example on syntetizing with Q-filter (forward Q-filtering) .

clear all

presets

global S4M

% Create a min-phase wavelet, zero-phase or mix-phase

wavelet=s_create_wavelet({’type’,’min-phase’},{’step’,1});

% Create constant-Q filters for Q = 150,100,80,60 and 1 sec travel time

qfilters=s_create_qfilter({’q’,[150,100,80,60]},{’times’,1000},{’step’,1});

% Convolve the constant-Q filters with the wavelet qwavelets=s_convolve(qfilters,wavelet);

% Prepend the original wavelet (since the original wavelet is shorter

% than the filtered ones the missing samples are, by default, set to

% NaN’s; this will give a warning in R 14, if the data set is plotted}

wavelets=s_append(wavelet,qwavelets);

% Give the data set “wavelets” a descriptive name wavelets.name=’Original wavelet and its Q-filtered versions’;

% Plot the result with the following conditions:

% plot only the first 200 ms;

% annotate traces 2 to 5 with the Q-value of the filter (the Q-value

% of the first trace is NaN as it is the original wavelet);

% fill the wiggle troughs with light gray

s_wplot(wavelets,{’times’,0,200},{’aindex’,2:5},{’annotation’,’Q’}, ...

{’trough_fill’,[0.7,0.7,0.7]});

% Display real and imaginary parts of the Fourier transform of

% the wavelets

ftwavelets=s_fft(wavelets,{’output’,’ft’});

s_compare(real(ftwavelets),imag(ftwavelets),{’times’,0,80})

stitle(‘Real (black) and imaginary (red) parts of the Fourier transform of the wavelets’)

%The package Seislab that I use is opensource. That means that the source code can be changed.

I will like to implement the Riccati-equation in the code and synthetize so transmission losses and multiples are introduced in the traces.(this will be done later.)

%Inversion will remove these effects and when a good inverse Q-filter is introduced, absorption is also removed.

%We see that lower Q gives more damping, as can be clearly seen on frequency plot. Both amplitude and frequency content is reduced considerably.

Application of inverse Q-filtering on test-data.

Our research must be tested on real data before we can see the usefulness of the methods. The Riccatiequation has been applied to real data in a very small extent up till now. So we introduce a test now.

Test data is given in Seislab on SEG format and this data will be deconvolved. Can easily syntethize a log for testing:

clear all

presets

global S4M

S4M.interactive=0; % Run without user intervention (turn-off some pop-up

% windows requesting user intervention)

% Create a synthetic log with a sonic made up from Gaussian noise;

depth=(5000:0.5:10000)';

nsamp=length(depth);

dt=150+25*randn(nsamp,1);

wlog=l_convert([depth,dt],{'depth','ft','Depth';'DTp','us/ft','P-sonic'});

wlog.name='Synthetic log';

% Add Kelley bushing elevation as a parameter

wlog=ds_add_parameter(wlog,84,{'KBE','ft','Kelly bushing elevation'});

% Add a trend to the sonic log

wlog=l_curve_math(wlog,'replace','dtp=dtp*(5000/depth)^0.333');

% Extract sonic curve from log

dtp=l_gc(wlog,'dtp');

% Use "dtp" to add a curve with P-velocity (curve mnemonic 'Vp')

wlog=l_curve(wlog,'add','Vp',1.0e6./dtp,'ft/s','P-velocity');

% Use an alternative method to add a density curve computed

% via Gardner's formula

wlog=l_curve_math(wlog,'add_ne','rho=0.23*vp^0.25','g/cm3','Density');

% Plot curves

l_plot(wlog) %#ok

%make reflection coefficient

%problem: how to make a reflection coefficient series from the wlog? Now I

%make it out of random noise

% Create 50 reflection coefficient series from random Gaussian noise

reflect=s_convert(randn(251,50),0,4);

% Create a Ricker wavelet

wavelet=s_create_wavelet({'type','Ricker'});

% Convolve reflection coefficient series with wavelet

synthetic=s_convolve(reflect,wavelet);

% Add a header (CDP) to the synthetic

synthetic=ds_header(synthetic,'add','CDP',101:150,'n/a','CDP number');

% Select the time range from 0 to 1000 ms

synthetic=s_select(synthetic,{'times',0,1000});

% Create a color plot of the synthetic

s_cplot(synthetic);

% Write the synthetic to a SEG-Y file

filename=fullfile(tempdir,'test.sgy');

write_segy_file(synthetic,filename);

% Read a subset of the SEG-Y file (include traces with even CDP's less than 121)

seismic=read_segy_file(filename,{'traces','CDP < 121 & mod(CDP,2) == 0'});

% Make a wiggle-trace plot of the seismic data set (specify peak fills

% and trough fills other than the default and annotate traces by CDP number)

s_wplot(seismic,{'trough_fill','gray'},{'peak_fill','red'},{'annotation','CDP'})

The wiggle-trace plot of the seismic data set shows fills and trough other than the default and annotate traces by CDP number.

% Apply an Ormsby band-pass filter with corner frequencies 0, 10,

% 20, 40 Hz to the seismic

fseismic=s_filter(seismic,{'ormsby',0,10,20,40});

% Plot the filtered seismic on top of the unfiltered seismic

s_compare(seismic,fseismic)

mytitle('Original (black) vs. filtered (red) seismic')

The bandpass-filtering in red shows some frequency-loss on corner frequency 0,10,20,40. We have plotted the spectra of the trace:

Before attempt to decon test data we must test the inverse Q-filter. Then we apply a forward Q-filter on some of the test data: