How to extract data from txt file and plot spectogram?
Afficher commentaires plus anciens
Hello,
I have a large sample of data (100000 samples) in txt file consisting of time and amplitude of the signal.
I want to plot spectogram. The code is not at all able to extract the .txt file at all along the path.
Kindly advice me as how to extract data from txt file when data is very large as it is in my case and hence to plot the spectogram.
The code is as below:
close all;
freq_sam = 100000; % samples in tt1 data
period_tt1 = 0.005; % time period
freq_tt1 = 1/period_tt1; % frequency
data = importdata('D:\Data_tokamak\data\obp.txt');
x_data = data(:,1);
y_data = data(:,2);
disp(x_data);
disp(y_data);
signal_tt1 = data;
fft_tt1 = fft(signal_tt1); % fourier transform
fft_tt1 = fftshift(fft_tt1); % performing fft shift
f = freq_sam/2*linspace(-1,1,freq_sam);
figure;
plot(f, fft_tt1);
Réponse acceptée
Star Strider
le 3 Juin 2023
I prefer the pspectrum function for these analyses. See specifically Spectrogram and Reassigned Spectrogram of Chirp for an example of returning the data and plotting it as well.
data = readmatrix('D:\Data_tokamak\data\obp.txt');
x_data = data(:,1);
y_data = data(:,2);
Ts = mean(diff(x_data))
Fs = 1/Ts;
figure
pspectrum(y_data, Fs, 'spectrogram') % Plot 'pspectrum' 'spectrogram'
[sp,fp,tp] = pspectrum(y_data,fs,"spectrogram"); % Return Values, Then Plot
figure
waterfall(fp,tp,sp')
set(gca,XDir="reverse",View=[60 60])
ylabel("Time (s)")
xlabel("Frequency (Hz)")
Be certain that the ‘x_data’ are sampled uniformly. The easiest way to do that is to calculate the standard deviation (std) of ‘diff(x_data)’ and compare it to ‘Ts’. It should be vanishingly small, and preferably less than 1E-6 of the ‘Ts’ value.
It would help to have ‘obp.txt’ to work with.
.
37 commentaires
Rahul
le 3 Juin 2023
I doubt if it is too large, because 1E+5 samples is reasonable. Without having the data to work with, I suspect the problem is the presence of Inf or NaN values in ‘y_data’.
The fillmissing function searches only for NaN (or NaT — ‘Not a Time’) values, so the first thing to do is to search for those values in ‘y_data’ and first replace them with NaN values then use fillmissing to interpolate them —
data = [(0:10).' [rand(4,1); Inf; rand(3,1); NaN; rand(2,1)]]
data(~isfinite(data)) = NaN
data = fillmissing(data, 'linear')
I would not remove the missing value rows because spectral analysis functions (with the exception only of the nufft function that I mention here only to note its presence) require (and assume) regular sampling, and removing the rows with missing values violates that assumption. (The interpolation method can be whatever you want, I chose 'linear' here for convenience.)
The slightly modified version of ‘data’ should then work.
.
Star Strider
le 4 Juin 2023
Modifié(e) : Star Strider
le 4 Juin 2023
The line:
data(~isfinite(data)) = NaN;
should not convert all the values to NaN. I demonstrated how it works earlier. In my code it does not need an if block to work correctly. It already detects all of them and then simply converts all the Inf values to NaN so that fillmissing can interpolate them.
Another option is:
data(isinf(data)) = NaN;
data = fillmissing(data, 'linear');
and then the rest of the code in my original Answer should work without problems.
Please share your data. Attach it using the ‘paperclip’ icon in the top toolbar. That is the easiest way to get this to work.
.
Rahul
le 4 Juin 2023
Thank you!
Try this —
filecontents = fileread('obp3n.txt') % Check File Format & Contents
data = readmatrix('obp3n.txt', 'HeaderLines',8) % Skip First 8 Header Lines
nonfinite_rows = nnz(~any(isfinite(data),2)) % Check To See If Any Rows Have Non-Finite Values
if nonfinite_rows > 0
data(~isfinite(data)) = NaN;
data = fillmissing(data, 'linear');
end
x_data = data(:,1);
y_data = data(:,2);
Ts = mean(diff(x_data))
Tsd = std(diff(x_data))
Fs = 1/Ts
[y_data,x_data] = resample(y_data, x_data, round(Fs)); % Resample To Constant Sampling Frequency
Fs = round(Fs);
figure
pspectrum(y_data, Fs, 'spectrogram') % Plot 'pspectrum' 'spectrogram'
colormap(turbo) % Introduced: R2020b
[sp,fp,tp] = pspectrum(y_data,Fs,"spectrogram"); % Return Values, Then Plot
figure
waterfall(fp,tp,sp')
set(gca,XDir="reverse",View=[60 60])
colormap(turbo) % Introduced: R2020b
ylabel("Time (s)")
xlabel("Frequency (Hz)")
figure
waterfall(fp,tp,sp')
set(gca,XDir="reverse",View=[60 60])
Ax = gca;
Ax.ZScale = 'log';
colormap(turbo) % Introduced: R2020b
ylabel("Time (s)")
xlabel("Frequency (Hz)")
It was necessary to change the readmatrix call arguments to import the file correctly. It now works, and gives reasonable results.
There were enough variation in the sampling intervals thaat I considered resampling the data to be worthwhile. That is not necessary, however I believe it improved the result.
.
Rahul
le 4 Juin 2023
Star Strider
le 4 Juin 2023
My pleasure!
You will probably need to use the fullfile function to create a file name that includes the relevant path information to your other files. Once you have the directory information, simply change the file name argument in fullfile to create the complete path, and then use that as the file name.
For this file here, that might go something like this:
f = @(filename) fullfile('https://www.mathworks.com', 'matlabcentral/answers/uploaded_files/1402564', filename);
file = f('obp3n.txt')
data = readmatrix(file, 'HeaderLines',8) % Skip First 8 Header Lines
That creates:
file = 'https://www.mathworks.com/matlabcentral/answers/uploaded_files/1402564/obp3n.txt'
You will have to experiment with it on your computer, although it would be easier if all your files were added to your MATLAB search path, then the fullfile function (and the anonymous funciton I created here to make its use more straightforward) would not be necessary. See: What Is the MATLAB Search Path? for those details.
If my Answer helped you solve your problem, please Accept it!
.
Rahul
le 5 Juin 2023
Star Strider
le 5 Juin 2023
No NaN values are generated.
These lines:
nonfinite_rows = nnz(~any(isfinite(data),2)) % Check To See If Any Rows Have Non-Finite Values
if nonfinite_rows > 0
data(~isfinite(data)) = NaN;
data = fillmissing(data, 'linear');
end
check to see if there are any non-finite values in the file, and if there are, it fills them so that the structure of the data remains intact. It does not linearize anything, it uses linear interpolation (other methods are possible) to fill any non-finite values. That avoids problems.
If you are certain that the data will never have any non-finite values, you can eliminate those lines. The code will run correctly without them.
.
Multiply ‘data(:,1)’ by 1E-3 and it works —
filecontents = fileread('obp3n.txt') % Check File Format & Contents
data = readmatrix('obp3n.txt', 'HeaderLines',8) % Skip First 8 Header Lines
data(:,1) = data(:,1) * 1E-3;
nonfinite_rows = nnz(~any(isfinite(data),2)) % Check To See If Any Rows Have Non-Finite Values
if nonfinite_rows > 0
data(~isfinite(data)) = NaN;
data = fillmissing(data, 'linear');
end
x_data = data(:,1);
y_data = data(:,2);
Ts = mean(diff(x_data));
Tsd = std(diff(x_data));
Fs = 1/Ts
[y_data,x_data] = resample(y_data, x_data, round(Fs)); % Resample To Constant Sampling Frequency
Fs = round(Fs);
figure
pspectrum(y_data, Fs, 'spectrogram') % Plot 'pspectrum' 'spectrogram'
colormap(turbo) % Introduced: R2020b
[sp,fp,tp] = pspectrum(y_data,Fs,"spectrogram"); % Return Values, Then Plot
figure
waterfall(fp,tp,sp')
set(gca,XDir="reverse",View=[60 60])
colormap(turbo) % Introduced: R2020b
ylabel("Time (s)")
xlabel("Frequency (Hz)")
figure
waterfall(fp,tp,sp')
set(gca,XDir="reverse",View=[60 60])
Ax = gca;
Ax.ZScale = 'log';
colormap(turbo) % Introduced: R2020b
ylabel("Time (s)")
xlabel("Frequency (Hz)")
Fixed!
.
Rahul
le 7 Juin 2023
I generally prefer to use elliptic filters. If you have R2018a or later, you can just do this:
y_data_filt = bandpass(y_data, [35 50], Fs, 'ImpulseResponse','iir');
Otherwise, this works —
data = readmatrix('obp3n.txt', 'HeaderLines',8) % Skip First 8 Header Lines
data(:,1) = data(:,1) * 1E-3; % Time In Milliseconds
x_data = data(:,1);
y_data = data(:,2);
Ts = mean(diff(x_data));
Tsd = std(diff(x_data));
Fs = 1/Ts
[y_data,x_data] = resample(y_data, x_data, round(Fs)); % Resample To Constant Sampling Frequency
Fs = round(Fs);
Fn = Fs/2;
Wp = [35 50]/Fn; % Passband Frequency (Normalised)
Ws = [0.98 1.02].*Wp; % Stopband Frequency (Normalised)
Rp = 1; % Passband Ripple
Rs = 60; % Passband Ripple (Attenuation)
[n,Wp] = ellipord(Wp,Ws,Rp,Rs); % Elliptic Order Calculation
[z,p,k] = ellip(n,Rp,Rs,Wp); % Elliptic Filter Design: Zero-Pole-Gain
[sos,g] = zp2sos(z,p,k); % Second-Order Section For Stability
figure
freqz(sos, 2^18, Fs) % Filter Bode Plot
set(subplot(2,1,1), 'XLim',[0 100]) % OPTTIONAL
set(subplot(2,1,2), 'XLim',[0 100]) % OPTTIONAL
y_data_filt = filtfilt(sos, g, y_data); % Filter Signal
figure
plot(x_data, y_data_filt)
grid
xlabel('Time')
ylabel('Amplitude')
figure
pspectrum(y_data_filt, Fs, 'spectrogram') % Plot 'pspectrum' 'spectrogram'
colormap(turbo) % Introduced: R2020b
[sp,fp,tp] = pspectrum(y_data,Fs,"spectrogram"); % Return Values, Then Plot
figure
waterfall(fp,tp,sp')
set(gca,XDir="reverse",View=[60 60])
colormap(turbo) % Introduced: R2020b
ylabel("Time (s)")
xlabel("Frequency (Hz)")
figure
waterfall(fp,tp,sp')
set(gca,XDir="reverse",View=[60 60])
Ax = gca;
Ax.ZScale = 'log';
colormap(turbo) % Introduced: R2020b
ylabel("Time (s)")
xlabel("Frequency (Hz)")
.
Rahul
le 8 Juin 2023
Star Strider
le 8 Juin 2023
That relates to the way the axes are scaled. They all go from 0 Hz (D-C) to 
Hz, equivalent to 
kHz. The first plot is scaled in kHz and goes from 0 to 
kHz, the others are scaled in Hz and go from 0 Hz to 
(equivalent to 
) Hz.
Hz, equivalent to kHz. The first plot is scaled in kHz and goes from 0 to kHz, the others are scaled in Hz and go from 0 Hz to (equivalent to ) Hz. So, they are all the same.
Rahul
le 11 Juin 2023
The freqz function with outputs will not also plot the Bode plot, so two separate calls to it are necessary here, however the second argument to it in this code limits the fft to only length 10 (2^3.3), while I usually use something on the order of 2^18.
The result is that will not be anywhere close to being an accurate Bode plot for this filter —
Fs = 2.0000e+05
Fn = Fs/2;
Wp = [35 50]/Fn; % Passband Frequency (Normalised)
Ws = [0.98 1.02].*Wp; % Stopband Frequency (Normalised)
Rp = 1; % Passband Ripple
Rs = 60; % Passband Ripple (Attenuation)
[n,Wp] = ellipord(Wp,Ws,Rp,Rs); % Elliptic Order Calculation
[z,p,k] = ellip(n,Rp,Rs,Wp); % Elliptic Filter Design: Zero-Pole-Gain
[sos,g] = zp2sos(z,p,k); % Second-Order Section For Stability
figure(4)
freqz(sos, 10,Fs) ; % Filter Bode Plot
% set(subplot(2,1,1), 'XLim',[0 100]) % OPTTIONAL
% set(subplot(2,1,2), 'XLim',[0 100]) % OPTTIONAL
[h,w]=freqz(sos, 10,Fs) ;
figure(5)
ax = gca;
ax.YLim=[0,100];
ax.XTick = 0:10:100;
plot(w, pi/180*unwrap(angle(h)));
xlabel('Frequency (Hz)');
ylabel('Phase (rad)');
I am not certain what you want to do here.
.
Rahul
le 11 Juin 2023
Star Strider
le 11 Juin 2023
Rahul
le 13 Juin 2023
Déplacé(e) : Star Strider
le 13 Juin 2023
Star Strider
le 13 Juin 2023
I have no idea how you generate ‘h’ (that previously was created by freqz) or how it relates to ‘Mydata’ that your code never uses or refers to after creating it. If you are supposed to filter it after reading it, you need to add a filtfilt call using the previously-created filter, and then do sometihing with the fitlered data.
There is no freqz call anywhere. If it is called before the loop, the code plots the same information 12 times.
I also changed the file name format descriptor to match the previous file names. MATLAB is case-sensitive, and the file names must match exactly.
I made some minor changes, however I cannot fix it because I do not know what you want to do with ‘Mydata’ or even what it is —
figure(5)
hold on
for i=1:12
myfilename = sprintf('obp%in.txt',i);
file = fullfile('D:\Data_tokamak\TT1 data',myfilename );
Mydata = readmatrix(file, 'HeaderLines',8);
ax = gca;
ax.YLim=[0,100];
ax.XTick = [0:10:100];
plot(w, pi/180*(angle(h)));
grid on;
xlabel('Frequency (Hz)');
ylabel('Phase (rad)');
end
hold off
This code honestly does not make sense to me.
.
Rahul
le 13 Juin 2023
Star Strider
le 13 Juin 2023
As always, my pleasure!
Rahul
le 13 Juin 2023
Star Strider
le 13 Juin 2023
Modifié(e) : Star Strider
le 13 Juin 2023
There are several problems with the first loop.
The first loop should be something like this —
f = @(filename)fullfile('D:\Data_tokamak\TT1 data',filename);
for k = 1:numfiles
myfilename = sprintf('obp%dn.txt', k);
mydata{k} = readmatrix(f(myfilename), 'HeaderLines',8);
mydata{k}(:,1) = mydata(:,1) * 1E-3;
end
That looks like it should work.
Using this construction, this needs to be changed to:
x_data(:,k) = mydata{k}(:,1);
y_data(:,k) = mydata{k}(:,2);
The rest of the code seems to be what has already been demonstrated to work, at least with the provided file.
I cannot test it beyond wiat I already wrote.
.
Rahul
le 14 Juin 2023
Star Strider
le 14 Juin 2023
‘Too many input arguments’
What line is throwing that error? (Please copy and paste all the red text from the Command Window.)
What arguments are you passing to that function call?
It is UCT-6 here (22:25) so I will deal with this in the morning.
.
Rahul
le 14 Juin 2023
Star Strider
le 14 Juin 2023
O.K.
mydata{k} = readmatrix(f(myfilename), 'HeaderLines',8);
That should work.
Rahul
le 14 Juin 2023
Star Strider
le 14 Juin 2023
Rahul
le 14 Juin 2023
Star Strider
le 14 Juin 2023
[y_datac{:,k},x_datac{:,k}] = resample(y_data(:,k), x_data(:,k), round(Fs));
That should work. The only complication is that you will have to use cell array indexing on ‘x_datac’ and ‘y_datac’ after that. It is not possible to store them in ‘x_data’ and ‘y_data’ because those were not originally cell arrays, and the sizes do not match, so they cannot be stored as numeric arrays.
.
The problem seems to be whatever ‘k’ is, since it works here for different values of ‘k’ —
Fs = 200000;
Fn = Fs/2;
Wp = [35 50]/Fn; % Passband Frequency (Normalised)
Ws = [0.98 1.02].*Wp; % Stopband Frequency (Normalised)
Rp = 1; % Passband Ripple
Rs = 60; % Passband Ripple (Attenuation)
[n,Wp] = ellipord(Wp,Ws,Rp,Rs); % Elliptic Order Calculation
[z,p,k] = ellip(n,Rp,Rs,Wp); % Elliptic Filter Design: Zero-Pole-Gain
[sos,g] = zp2sos(z,p,k); % Second-Order Section For Stability
k = 5;
y_datac{:,k} = randn(50,1);
y_data_filt = filtfilt(sos, g, y_datac{:,k});
CheckSize = size(y_data_filt)
So ‘k’ must be an integer greater than zero. (The code does not use logical indexing, so we can rule that out.)
.
Rahul
le 15 Juin 2023
Star Strider
le 15 Juin 2023
Are you certain this is correct:
myfilename = sprintf('OBP%dN.txt', k);
MATLAB is case-sensitive, and the file names must match exactly, even with respect to case, and the other file names were all lower-case. If that is the problem, the code should throw an error about not being able to find the file. Are there any errors?
So the first thing to do is to check to see what ‘mydata’ contains. Be certain that it is not empty.
Plus de réponses (1)
Diwakar Diwakar
le 3 Juin 2023
Try this code. May be this code will help you.
% Specify the path to your text file
file_path = 'path/to/your/file.txt';
% Read the data from the text file
data = dlmread(file_path);
% Extract the time and amplitude columns
time = data(:, 1);
amplitude = data(:, 2);
% Set the parameters for the spectrogram
window_size = 256; % Size of the window for the spectrogram
overlap = 128; % Overlap between consecutive windows
% Compute the spectrogram
[~, frequency, time_spectrogram, power] = spectrogram(amplitude, window_size, overlap, [], 'onesided', 'yaxis');
% Plot the spectrogram
figure;
imagesc(time_spectrogram, frequency, 10*log10(power));
axis xy; % Flip the y-axis direction
colormap jet; % Choose a colormap
colorbar; % Add a colorbar for reference
xlabel('Time');
ylabel('Frequency');
title('Spectrogram');
4 commentaires
Rahul
le 3 Juin 2023
Diwakar Diwakar
le 4 Juin 2023
data = readmatrix(file_path);
If you are using MATLAB R2019b or later, you can use the readmatrix function instead of dlmread
Rahul
le 4 Juin 2023
Star Strider
le 4 Juin 2023
@Rahul — We already covered this and I told you how to fix it.
If we are going to solve it, you need to attach your data.
Catégories
En savoir plus sur Multirate Signal Processing dans Centre d'aide et File Exchange
le 3 Juin 2023
le 15 Juin 2023
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
