How can I plot some fft data in a different way?

9 views (last 30 days)
Dear all,
I usual plot spectrograms plotting time vs. frequency vs. amplitude. I have collected some records through a microphone and each of the recordings are named usign a parameter, so called 'ϕ'. I would like to plot a single spectrogram similar to the one in the attached picture containing my data. File structure is 2 colums (time and amplitude) and 10000 rows.
Could you please help me making this "ϕ vs. frequency vs. amplitude" plot by using the two signal files attached?
I would sincerly appreciate it since I am litterally freaking out trying to fixing this issue.
  4 Comments
Yazan
Yazan on 22 Jul 2021
Edited: Yazan on 22 Jul 2021
Where are the values of ϕ? In the txt files, I can see time and amplitude. Also, what do you mean by close to each others? on the same axes (appending the first signal with the second?)? Or two subplots?
Francesco Pignatelli
Francesco Pignatelli on 22 Jul 2021
Edited: Francesco Pignatelli on 22 Jul 2021
Values of ϕ are in the file name: 0.66 and 0.68,respectively.
When I say close to each other I mean attached/stitched. I attach an example regarding what I would like to have.

Sign in to comment.

Accepted Answer

Yazan
Yazan on 22 Jul 2021
Edited: Yazan on 22 Jul 2021
The first column in the txt files you provided is a DateTime vector, so again, I am not able to understand how ϕ is defined. I will assume that the seconds in the DateTime vector represent ϕ.
clc, clear, close all
% read data
opts = detectImportOptions('CH4_NO_PILOT_RE_10K_PHI_066.txt');
opts.VariableNames = {'Time', 'Amp'};
opts.DataLines = [5 inf];
ch1 = readtable('CH4_NO_PILOT_RE_10K_PHI_066.txt', opts);
ch2 = readtable('CH4_NO_PILOT_RE_10K_PHI_068.txt', opts);
% get the time instants at which the signal is sampled
ch1.Time = seconds(ch1.Time - ch1.Time(1));
ch2.Time = seconds(ch2.Time - ch2.Time(1));
% sampling period
Ts = ch1.Time(2) - ch1.Time(1);
% downsample by a factor of 5 to reduce the number of datapoints
% this part can be skipped
sig1 = resample(ch1.Amp, 1, 5);
sig2 = resample(ch2.Amp, 1, 5);
Ts = Ts*5;
% histogram computation
[sp1, f, t] = pspectrum(sig1, 1/Ts, 'spectrogram', 'OverlapPercent', 25, 'Leakage', 0.8);
sp2 = pspectrum(sig2, 1/Ts, 'spectrogram', 'OverlapPercent', 25, 'Leakage', 0.8);
% append the second histogram with the first
spTot = pow2db([sp1, sp2]); t = [t; t]; Nt = length(t)/2;
% plot
figure,
ax = axes(gcf);
imagesc(1:2*Nt, f./1000, spTot);
cb = colorbar; axis xy
% label the x-axis as phi (usually it should be time)
xlabel('\phi', 'FontSize', 15),
% label the other axes
ylabel('Frequency (KHz)', 'FontSize', 15)
cb.Label.String = 'Power (dB)';
% change the x-axis ticks
ax.XTick = Nt*(0.5:0.5:2);
ax.XTickLabel = arrayfun(@(x) num2str(x), t(ax.XTick), 'UniformOutput', false);
% vertical line to separate the two histograms
xline(Nt, 'LineWidth', 3, 'LineStyle', '--', 'Color', 'r')
% annotate the figure
annotation('textbox', 'Position', [0.18,0.8,0.22,0.1], 'String', 'Histogram 1',...
'FontSize', 15)
annotation('textbox', 'Position', [0.5,0.8,0.22,0.1], 'String', 'Histogram 2',...
'FontSize', 15)
  3 Comments
Yazan
Yazan on 22 Jul 2021
Use the following
% change the x-axis ticks
ax.XTick = Nt*([0.5 1.5]);
ax.XTickLabel = [{'\phi = 0.66'}, {'\phi = 0.68'}];
Francesco Pignatelli
Francesco Pignatelli on 22 Jul 2021
Now it is perfect! Thank you very much for your help Yazan, I really appreciate it :)
I hope to "meet" you again in this MATLAB comunity. Bye :)

Sign in to comment.

More Answers (1)

Chunru
Chunru on 22 Jul 2021
% Assume there are 3 spec (corresponding to phi)
nf = 20; nt=15; nspec = 3;
sp = randn(nf, nt, nspec);
% reshape the spectrogram to 2D
sp = reshape(sp, [nf, nt*nspec]);
% multiple spectrograms
imagesc(sp);
% annotation
h = gca;
h.XTick = 0.5+(0:nspec-1)*nt;
h.YTick = [];
grid on
phi = [1.2 3.4 5.6];
h.XTickLabel = string(phi);
xlabel('\phi')
  1 Comment
Francesco Pignatelli
Francesco Pignatelli on 22 Jul 2021
Hello Chunru, thank you very much for your efforts, I really appreciate it :)

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!