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

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.

[Yazan]

Define for us first the relation between time and ϕ. As you said, the spectrogram helps you plot the instantaneous spectra. To plot the instantaneous spectra as a function of ϕ, you need to define the relation between ϕ and time.

[Francesco Pignatelli]

Well, there's not a mathematical or physical relation between time and ϕ. The latter is just a parameter useful for me to recognize the differences between the two cases.

Probably I didn't explain myself in the best way: I would like to plot the spectrogram of my two signal close to each other but I would like to plot values of ϕalong the x-axis (similarly to the attached picture).

[Yazan]

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]

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.

Answer 1

[Yazan]

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)

[Francesco Pignatelli]

I have tested your code and it seems working regularly out of the following two lines:

[sp1, f, t] = pspectrum(sig1, 1/Ts, 'spectrogram', 'OverlapPercent', 25, 'Leakage', 0.8);

sp2 = pspectrum(sig2, 1/Ts, 'spectrogram', 'OverlapPercent', 25, 'Leakage', 0.8);

since my Matlab complains with the following error message:

"The denominator in matrix division for duration arrays must be a real scalar double value."

that I have fixed just by not using Ts as duration but as a value.

The latest thing regarding ϕ: what I want is something similar to the following:

So that it is possible to distinguish the first dataset from the second one.

[Yazan]

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]

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 :)

Answer 2

[Chunru]

% 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')

[Francesco Pignatelli]

Hello Chunru, thank you very much for your efforts, I really appreciate it :)