I need to draw a spectrogram for given time period.

Hi, I need to draw a spectrogram of two signal. I have seen spectrogram function in help but could not relate it with my problem. Can anyone help me in function how can I get the spectrogram for specific time periods. Thanks.

clear all

load CO2data % Historical CO2 data from Harris

ageCO2=CO2data.years;

fCO2=CO2data.fCO2;

timeCO2=1950-ageCO2;

size(timeCO2)

timeCO2=flipud(timeCO2);

fCO2=flipud(fCO2);

timeCO2(565)=[]; % remove data 565 because it is double and gives an error in interp1

fCO2(565)=[]; % remove data 565 because it is double

figure(101)

plot(timeCO2/1000,fCO2)

xlabel('time (1000 years)')

ylabel('CO2 level (ppm)')

title('Historical CO2 level in the atmosphere (ppm)')

load Tdata % Historical temp data from Harris

ageT=Tdata.years;

Temp=Tdata.T;

timeT=1950-ageT;

size(timeT)

timeT=flipud(timeT);

Temp=flipud(Temp);

figure(102)

plot(timeT/1000,Temp)

xlabel('time (1000 years)')

ylabel('Temperature anomaly (C)')

title('Historical temperature anomaly (C)')

% Interpolation

time_int=-10000:1000; % interpolated time (want to draw spectrogram for this specific time period)

fTemp_int=interp1(timeT,Temp,time_int,'linear'); % interpolated Temp

fCO2_int=interp1(timeCO2,fCO2,time_int,'linear'); % interpolated CO2

figure(101)

plot(timeCO2/1000,fCO2,'b-')

hold on

plot(time_int/1000,fCO2_int,'r--')

hold off

xlabel('time (1000 years)')

ylabel('CO2 level (ppm)')

title('Historical CO2 level in the atmosphere (ppm)')

legend('original CO2 data','interpolated CO2 data')

figure(102)

plot(timeT/1000,Temp,'b-')

hold on

plot(time_int/1000,fTemp_int,'r--')

hold off

xlabel('time (1000 years)')

ylabel('Temperature anomaly (C)')

title('Historical temperature anomaly (C)')

legend('original temp data','interpolated temp data')

% Compute cross correlation

crosscorrelation=xcorr(fCO2_int-mean(fCO2_int),fTemp_int-mean(fTemp_int));

M=length(time_int);

figure(103)

plot(-M+1:M-1,crosscorrelation)

xlabel('time lag in years')

title('cross correlation function')

[Mathieu NOE]

hello

seems the topic of your code is about crosscorrelation, what does spectrogram have to do here ? this will give you the frequency / amplitude distribution vs time for each individual signal but no information about their possible inter correlation ...

[TEJASHKUMAR PATEL]

Hello,

I agree what you have said. In my case it is about the absolute time vs time lag in year. So I do not know much but by spectrogram I think I can get the time lag where the cross correlation is maximum from the data of past 800,000 years of temperature and CO2. I am not sure it will help or not but I thought so so I need to draw spectrogram. I tried to do it by imagesc but I do not know exactly how can I implement it.

[Mathieu NOE]

This is an example of spectrogram, either with the native matlab function, or my own implementation, supplied inside the code

fs=5000; % Sampling rate
N=1000; % Number of data points
t=[0:1:N-1]/fs; 
y=sin(500*pi*t)...
    +2*sin(500*1.5*pi*t)...
    +3*sin(500*3*pi*t)...
    +4*sin(500*4*pi*t);
 
figure(1); 
windowsize = 200; 
window = boxcar(windowsize); 
nfft = windowsize; 
noverlap = windowsize-1; 
% [S,F,T] = spectrogram(y,window,noverlap,nfft,fs); % S = SPECTROGRAM(X,WINDOW,NOVERLAP,NFFT,Fs)
[S,F,T] = myspecgram(y, fs, nfft, 0.75); % overlap is 75% of nfft here

imagesc(T,F,(abs(S))) 
colorbar('vert');
set(gca,'YDir','Normal') 
xlabel('Time (secs)') 
ylabel('Freq (Hz)') 
title('Short-time Fourier Transform spectrum')
colormap('jet');

function  [fft_specgram,freq_vector,time] = myspecgram(signal, Fs, nfft, Overlap)
% FFT peak spectrogram of signal  (example sinus amplitude 1   = 0 dB after fft).
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer overlap % (between 0 and 0.95)
signal = signal(:);
samples = length(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,1);
    s_tmp((1:samples),:) = signal;
    signal = s_tmp;
    samples = nfft;
end

% window : hanning
window = hanning(nfft);
window = window(:);

%    compute fft with overlap 
 offset = fix((1-Overlap)*nfft);
 spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
%     % for info is equivalent to : 
%     noverlap = Overlap*nfft;
%     spectnum = fix((samples-noverlap)/(nfft-noverlap));	% Number of windows

    % main loop
    fft_specgram = [];
    for ci=1:spectnum
        start = (ci-1)*offset;
        sw = signal((1+start):(start+nfft)).*window;
        fft_specgram = [fft_specgram abs(fft(sw))*4/nfft];     % X=fft(x.*hanning(N))*4/N; % hanning only 
    end

% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
fft_specgram = fft_specgram(select,:);
freq_vector = (select - 1)*Fs/nfft;

% time vector 
% time stamps are defined in the middle of the buffer
time = ((0:spectnum-1)*offset + round(offset/2))/Fs;

end

[TEJASHKUMAR PATEL]

Thanks for your time sir. I will check and will try to implement in my code.