How to find wavelength or wavenumber from signal? (MATLAB)

Hello

I am working on my master thesis and it is connected with signal processing.

I have a numerical model of plate in which i have got an extortion as a 5 periodic sin wave.

I have collected displacements in Y direction from all points on surface of my plate (400 points) in one time increment. I put it into Matlab as a vector:

data = load('jednachw100.txt');  
data=reshape(data,[],1);
y=fft2(data);
plot(abs(fft(data)));

and i have recived something like this:

In x axis i have my all sampling points which is 400(all nodes from surface of plate). My task is to get a wavenumber in x-axis but i can't do this.

Thaks for help in advance.

[Image Analyst]

What's your definition of wavenumber?

[Michal G]

definition of wavenumber is k=2pi/wavelength

Answer 1

[Walter Roberson]

Untested

numpoint = length(data);
curtick = get(gca, 'XTick');
lefthalf = curtick(curtick <= numpoint/2);
leftwaveL = numpoint ./ (lefthalf - 1);
righthalf = curtick(curtick > numpoint/2);
rightwaveL = numpoint ./ (numpoint - righthalf + 1);
waveL = [leftwaveL(:); rightwaveL(:)];
ticklab = cellstr(str2num(waveL));
set(gca, 'xticklabel', ticklab);

[Michal G]

Thanks for reply Walter but

Something is wrong...

??? Error using ==> cellstr at 34 Input must be a string.

Error in ==> test_fft at 17 ticklab = cellstr(str2double(waveL));

[Walter Roberson]

ticklab = cellstr(num2str(waveL));

[Michal G]

I think i should try something connected with this (n = 2^nextpow2(L);)

http://www.mathworks.com/help/matlab/ref/fft.html

[Walter Roberson]

I missed that data was the original data not the fft result. And you have no need of fft2 there.

datafft = fft(data);
plot(abs(datafft));
numpoint = length(data);
curtick = get(gca, 'XTick');
lefthalf = curtick(curtick <= numpoint/2);
leftwaveL = numpoint ./ (lefthalf - 1);
leftwaveL(lefthalf == 0) = inf;
righthalf = curtick(curtick > numpoint/2);
rightwaveL = numpoint ./ (numpoint - righthalf + 1);
waveL = [leftwaveL(:); rightwaveL(:)];
ticklab = cellstr(num2str(waveL));
set(gca, 'xticklabel', ticklab);

This is based upon the current tick positions, so that the wavelengths are shown at equal intervals. If you prefer "nice" but irregularly spaced numbers, then:

L = [1 2 3 4 8 16];
datafft = fft(data);
plot(abs(datafft));
numpoint = length(data);
leftpos = (numpoint-1)./(2.*L);
leftlab = L;
rightpos = numpoint - leftpos(2:end);
rightlab = L(2:end);
pos = [flipud(leftpos(:)); rightpos(:)];
labs = [flipud(leftlab(:)); rightlab(:)];
set(gca, 'xtick', pos);
set(gca, 'xticklabel', labs);

[Michal G]

So i am supposed to see the wavelength on x-axis? Or maybe wavenumber?

[Michal G]

I can add that i measured the wavelength on my numerical model and i know that the wavelength should be 0.0085 m.

[Walter Roberson]

The labeling I do here is samples per cycle, the wavelength. To convert that to a wavelength in meters we would have to know the total width.

Is the distance between your samples 0.425 mm, total length 17 cm? I estimate the peak is about 20 bins in, call it 21, which would be 20 cycles over the distance D (total length of the sampling area). If it is 20 cycles is 8.5 mm then your points are about 0.425 mm apart.

So given D, the total length in meters, the points I label with L in the code should have a wavelength of L * D/length(data).

My eyes keep closing here so before I change the code to label in those terms, I will wait for you to correct my estimates.

[Michal G]

Total width of my plate is 200mm = 0.2 m. Distance between samples is 0.5mm. And i have got 401 nodes exactly.

[Michal G]

What i hev measured on the model is about 8.5 mm not exactly 8.5mm.

[Walter Roberson]

Which is the constant, the 0.2 m, or the 0.5 mm? With 401 points, 0.5 mm would lead to 0.2005 m. If the 0.2 is correct then 401 points could be about 0.499 mm apart and if the peak is 17 cycles/distance that would be about 8.48 mm for the peak.

[Michal G]

plate is 200 mm. I have 400 elements, that means that i have 401 nodes. So the right value between points is 0.499mm. The pick is 23 cycles from left side. And constatn is of course 0.2 m.

[Michal G]

If u want i could send u my txt with my measurements.

[Walter Roberson]

Ummm. With

datafft = fft(data);

could you tell me the output of idx for

[~, idx] = max(abs(datafft(2:200)));

[Michal G]

22

[Walter Roberson]

You can attach the .txt file here, it isn't long. If you need to you can zip it and attach the .zip

[Walter Roberson]

It is 5 am here. I am going to bed now.

[Michal G]

ok, i have attached txt file. Just to remember. Txt consist y-dir displacements of all points from surface of a plate in one time increment.

[Michal G]

I think Walter that we should use k=2*pi/wavelength. That is my definition of wavenumber

[Walter Roberson]

I do not see the txt file attached. After you select a file, you have to click to Attach File.

[Michal G]

I can't add any more files today:( i sent yopu a private message with my vector.

Answer 2

[mm99]

Hi!

I am doing similar thing, but have some trouble with sine extortion. I get an error: "Unable to perform assignment because the size of the left side is 1-by-1 and the size of the right side is 1-by-101." This 1-by-1 is the point in which I want the sine to start and 1-by-101 is this sine with hann window. Do you maybe know what is wrong ?

Thanks in advance!