How to extract data from txt file and plot spectogram?

29 views (last 30 days)
Rahul
Rahul on 3 Jun 2023
Commented: Star Strider on 15 Jun 2023
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);

Sign in to answer this question.

Accepted Answer

Star Strider
Star Strider on 3 Jun 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 Comments
Show 35 older comments
Rahul
Rahul on 15 Jun 2023
Hi,
The value of k in my code is never initialized to zero, its going from 1 to numfiles(=12).
I did as you adviced. But none of the plots are getting executed.
Even the spectrogram is also not getting executed.
kindly advice.
with regards.
% Phase (in radians) versus Frequency (in Hz) plot for
% multiple input data from txt files
numfiles = 12;
x_datac{:,1}=[];x_datac{:,2}=[];x_datac{:,3}=[];x_datac{:,4}=[];x_datac{:,5}=[];x_datac{:,6}=[];x_datac{:,7}=[];
x_datac{:,8}=[];x_datac{:,9}=[];x_datac{:,10}=[];x_datac{:,11}=[];x_datac{:,12}=[];
y_datac{:,1}=[];y_datac{:,2}=[];y_datac{:,3}=[];y_datac{:,4}=[];y_datac{:,5}=[];y_datac{:,6}=[];y_datac{:,7}=[];
y_datac{:,8}=[];y_datac{:,9}=[];y_datac{:,10}=[];y_datac{:,11}=[];y_datac{:,12}=[];
mydata = cell(1, numfiles);
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{k}(:,1) * 1E-3;
end
for k=1:numfiles
nonfinite_rows{k} = nnz(~any(isfinite(mydata{k}),2)); % Check To See If Any Rows Have Non-Finite Values
if nonfinite_rows{k} > 0
data(~isfinite(mydata{k})) = NaN;
data = fillmissing(mydata{k}, 'linear');
end
x_data(:,k) = mydata{k}(:,1);
y_data(:,k) = mydata{k}(:,2);
Ts = mean(diff(x_data(:,k)));
Tsd = std(diff(x_data(:,k)));
Fs = 1/Ts;
[y_datac{:,k},x_datac{:,k}] = resample(y_data(:,k), x_data(:,k), round(Fs)); % Resample To Constant Sampling Frequency and using cell arrays
end
Fs = round(Fs);
Fn = Fs/2;
figure(1)
for k=1:numfiles
pspectrum(y_datac{:,k}, Fs, 'spectrogram'); % Plot 'pspectrum' 'spectrogram'
colormap(turbo); % Introduced: R2020b
[sp,fp,tp] = pspectrum(y_datac{:,k},Fs,"spectrogram"); % Return Values, Then Plot
end
wp = [39.8 40.2]/Fn; % Passband frequency Normalized
ws = [0.98 1.02].*wp; % Stopband frequency Normalized
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 selection for stability
figure(2)
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
[h w] = freqz(sos, 10,Fs) ;
%Implementation of filtfilt function without plotting
k=12;
y_data_filt = filtfilt(sos, g, y_datac{:,k}); % Filter Signal
% Plot for phase (in radians) versus frequency (in hz)
for k=1:numfiles
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)');
hold on
end
Star Strider
Star Strider on 15 Jun 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.

Sign in to comment.

More Answers (1)

Diwakar Diwakar
Diwakar Diwakar on 3 Jun 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 Comments
Show 2 older comments
Rahul
Rahul on 4 Jun 2023
Its not working
it gives error:
Error using pwelch
Expected x to be finite.
Plz advice.
rgds
Star Strider
Star Strider on 4 Jun 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.

Sign in to comment.

Categories

Find more on Multirate Signal Processing in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!