Extract signal values from a data file

Hello all,

I am trying to extract some basic signal values from a data file that I have which is a plot of a sinosodial wave vs length:

data = dlmread('data_30.txt');

z = data(:,1);

y = data(:,2);

figure(1)

plot(z,y,'b')

xlabel('Spanwise length')

ylabel('Velocity')

I want to identify and match the following:

% Fs=800;

% Tf=2;

% t=0:1/Fs:Tf;

% f=[40 75];

% Amp=[4.5 9.22];

% sigma=1.33;

% y=Amp(1)*exp(j*2*pi*t*f(1))+Amp(2)*exp(j*2*pi*t*f(2));

% N=(sigma/sqrt(2))*(randn(size(t))+j*randn(size(t)));

% y=y+N;

% fy=FFT(y,Fs);

The above code is part of the code:

PSD (Power Spectral Density), and Amplitude Spectrum with adjusted FFT

Version 1.3.0.0 (63.8 KB) by Youssef Khmou

I know that the length in my case is z, Can you guide me on how to identify the above values (Fs, Tf, f, Amp, N, etc..)?

The data file is attacehd, I appreciate your input.

Thanks!

Answer 1

[Star Strider]

I am not certain what ‘identify and match’ means. It is certainly possible to filter the 40 and 75 Hz frequencies individually, and their both defining a bandwidth, with a bandpass filter (introduced in R2018a, so you should have it) —

data = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1236212/data_30.txt');

z = data(:,1);
y = data(:,2);

Fs = 1/mean(diff(z))

Fs = 500

f=[40 75];

figure(1)
plot(z,y,'b')
xlabel('Spanwise length')
ylabel('Velocity')

Hz40 = bandpass(y, [-1 1]+f(1), Fs, 'ImpulseResponse','iir');
Hz75 = bandpass(y, [-1 1]+f(2), Fs, 'ImpulseResponse','iir');
Hz4075 = bandpass(y, f, Fs, 'ImpulseResponse','iir');

figure
subplot(3,1,1)
plot(z, Hz40)
grid
ylabel('Amplitude')
title('40 Hz')
subplot(3,1,2)
plot(z, Hz75)
grid
ylabel('Amplitude')
title('75 Hz')
subplot(3,1,3)
plot(z, Hz4075)
grid
xlabel('z')
ylabel('Amplitude')
title('40-75 Hz')

Experiment to get different results.

.

[Hussein Kokash]

First, thank you Star Strider for your valuable input, what I meant is using data similar to what I attached, how can I calculate the Power Spectral Density, speciialy that they are already calculated values and the code that I was looking at (mentioned earlier), deals with what you input such as:

% Fs=800;

% Tf=2;

% t=0:1/Fs:Tf;

% f=[40 75];

% Amp=[4.5 9.22];

% sigma=1.33;

% y=Amp(1)*exp(j*2*pi*t*f(1))+Amp(2)*exp(j*2*pi*t*f(2));

% N=(sigma/sqrt(2))*(randn(size(t))+j*randn(size(t)));

% y=y+N;

% fy=FFT(y,Fs);

Thanks again!

[Star Strider]

My pleasure!

I generally use the pwelch function to calculate the power spectral density (PSD). I have no idea what the code snippet you quoted does, so I have no idea how to incorporate it into my code or this code. If you want to filter it before computing the PSD, that is certainly an option.

You can also get the appropriate values by indexing into them in the output of pwelch. I do that here using ismembertol to get the indices and then present the results in two tables, the first with the absolute values and the second with decibels —

data = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1236212/data_30.txt');
z = data(:,1);
y = data(:,2);
Fs = 1/mean(diff(z))

Fs = 500

figure
pwelch(y)

[pxx,f,pxxc] = pwelch(y,[],[],[],Fs,'ConfidenceLevel',0.95);

fidx = ismembertol(f, [40 75], 0.01);
freqs = f(fidx);
PSDs = pxx(fidx);
PSDci = pxxc(fidx,:);

PSDTable = table(freqs,PSDs,PSDci)

PSDTable = 4×3 table
    freqs        PSDs                PSDci          
    ______    __________    ________________________

    39.062    1.2498e-10    6.9326e-11     2.895e-10
    41.016    1.3974e-10    7.7511e-11    3.2368e-10
    74.219    2.7447e-11    1.5224e-11    6.3575e-11
    76.172    1.9644e-11    1.0896e-11      4.55e-11

PSDdBTable = table(freqs,pow2db(PSDs),pow2db(PSDci),'VariableNames',{'freqs','PSDs (dB)','PSDci (dB)'})

PSDdBTable = 4×3 table
    freqs     PSDs (dB)        PSDci (dB)    
    ______    _________    __________________

    39.062     -99.031     -101.59    -95.384
    41.016     -98.547     -101.11    -94.899
    74.219     -105.62     -108.17    -101.97
    76.172     -107.07     -109.63    -103.42

figure
plot(f, pow2db(pxx), 'DisplayName','PSD')
hold on
plot(f, pow2db(pxxc), '--r', 'DisplayName','95% CI')
hold off
grid
xlabel('Frequency (Hz)')
ylabel('PSD (dB/Hz)')
legend('Location','best')

Experiment to get different results.

.

[Hussein Kokash]

I will test it accordingly, very infromative.

Thank you! Stay well!

[Star Strider]

As always, my pleasure!

[Hussein Kokash]

Hello , sorry for the trouble,

I have used this code but extended it so that it covers a 2D plane of coordinates (x, y) with their velocities (U, V)

clear all; clc;%close all;

path = pwd; % put current path here

S = dir(fullfile(path,'data_*.txt'));

addpath('/media/hf9098/easystore/MATLAB/POD-umich')

R = dir(fullfile(path,'data_*.txt'));

S = natsortfiles(R);

if iscell(S) == 0;

S = (S);

end

for k = 1:numel(S)

folder = S(k).folder;

filename = S(k).name;

F = fullfile(S(k).folder,S(k).name);

data = dlmread(F);

x = data(:,1); % x coordinate

y = data(:,2); % y coordinate

U(:,k) = data(:,4); % u velocity

V(:,k) = data(:,5); % v velocity

Fs(k) = 1/mean(diff(x));

figure(1)

pwelch(V)

[pxx,f,pxxc] = pwelch(V,[],[],[],Fs,'ConfidenceLevel',0.95);

fidx = ismembertol(f, [40 75], 0.01);

freqs = f(fidx);

PSDs = pxx(fidx);

PSDci = pxxc(fidx,:);

PSDTable = table(freqs,PSDs,PSDci)

PSDdBTable = table(freqs,pow2db(PSDs),pow2db(PSDci),'VariableNames',{'freqs','PSDs (dB)','PSDci (dB)'})

figure(2)

plot(f, pow2db(pxx), 'DisplayName','PSD')

hold on

plot(f, pow2db(pxxc), '--r', 'DisplayName','95% CI')

hold off

grid

xlabel('Frequency (Hz)')

ylabel('PSD (dB/Hz)')

legend('Location','best')

end

I am getting this error:

Error using table (line 307)

'PSDs (dB)' is not a valid variable name.

Error in testpwelch (line 36)

PSDdBTable = table(freqs,pow2db(PSDs),pow2db(PSDci),'VariableNames',{'freqs','PSDs (dB)','PSDci (dB)'})

Any thoughts?

Thanks!

[Star Strider]

As always, my pleasure!

The ability to use table variables that were not valid MATLAB variable names was introduced after R2018b. I forgot about that. Use: 'PSDs_dB' and such instead (or create your own variable names, since there is nothing particularly sacred about the ones I chose). Anything that is a valid MATLAB variable name should work.

Try this instead —

PSDdBTable = table(freqs,pow2db(PSDs),pow2db(PSDci),'VariableNames',{'freqs','PSDs_dB','PSDci_dB'})

That should run without error.

.