I am an undergraduate and new to matlab. i am doing my project in audio signal processing. can someone help me in drawing two frequency – time graph and compare them to mark the dissimilarities……

Thank you in advance…

Thenuja.J

Skip to content
# MATLAB: How to compare two frequency time graphs

#### Related Solutions

comparinggraphs

I am an undergraduate and new to matlab. i am doing my project in audio signal processing. can someone help me in drawing two frequency – time graph and compare them to mark the dissimilarities……

Thank you in advance…

Thenuja.J

hello

this code will do fft / spectrogram and bandpass filtering - adapt the corner frequencies to your own needs !

`clcclear all%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%`

% FFT parameters

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%NFFT = 4096; %

OVERLAP = 0.5;% spectrogram dB scale

spectrogram_dB_scale = 80; % dB range scale (means , the lowest displayed level is XX dB below the max level)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% options

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if you are dealing with acoustics, you may wish to have A weighted

% spectrums

% option_w = 0 : linear spectrum (no weighting dB (L) )

% option_w = 1 : A weighted spectrum (dB (A) )

option_w = 1;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% load signal

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% data

% [data,Fs] = audioread('myWAVaudiofile.wav');

[data,Fs] = audioread('Immigran_Hiss.wav');channel = 1;signal = data(:,channel);samples = length(signal);%% decimate (if needed)

% NB : decim = 1 will do nothing (output = input)

decim = 1;if decim>1 signal = decimate(signal,decim); Fs = Fs/decim;end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% display 1 : averaged FFT spectrum

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[freq, sensor_spectrum] = myfft_peak(signal,Fs,NFFT,OVERLAP);% convert to dB scale (ref = 1)

sensor_spectrum_dB = 20*log10(sensor_spectrum);% apply A weigthing if needed

if option_w == 1 pondA_dB = pondA_function(freq); sensor_spectrum_dB = sensor_spectrum_dB+pondA_dB; my_ylabel = ('Amplitude (dB (A))');else my_ylabel = ('Amplitude (dB (L))');endfigure(1),plot(freq,sensor_spectrum_dB,'b');gridtitle(['Averaged FFT Spectrum / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);xlabel('Frequency (Hz)');ylabel(my_ylabel);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% display 2 : time / frequency analysis : spectrogram demo

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[sg,fsg,tsg] = specgram(signal,NFFT,Fs,hanning(NFFT),floor(NFFT*OVERLAP)); % FFT normalisation and conversion amplitude from linear to dB (peak)

sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg)); % NB : X=fft(x.*hanning(N))*4/N; % hanning only

% apply A weigthing if neededif option_w == 1 pondA_dB = pondA_function(fsg); sg_dBpeak = sg_dBpeak+(pondA_dB*ones(1,size(sg_dBpeak,2))); my_title = ('Spectrogram (dB (A))');else my_title = ('Spectrogram (dB (L))');end% saturation of the dB range :

% saturation_dB = 60; % dB range scale (means , the lowest displayed level is XX dB below the max level)

min_disp_dB = round(max(max(sg_dBpeak))) - spectrogram_dB_scale;sg_dBpeak(sg_dBpeak<min_disp_dB) = min_disp_dB;% plots spectrogram

figure(2);imagesc(tsg,fsg,sg_dBpeak);colormap('jet');axis('xy');colorbar('vert');gridtitle([my_title ' / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(fsg(2)-fsg(1)) ' Hz ']);xlabel('Time (s)');ylabel('Frequency (Hz)');%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% removal of hiss between t = 0 and t = 1.68 s

% use of bandpass filter to remove noise below 500 and above 2500 Hz

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%data_filtered = data;t_min = 0;t_max = 1.68;ind = 1+fix(t_min*Fs:t_max*Fs-1);% keep only signal content extracted by bandpass filters

f_low = 500;f_high = 2500;N = 4;% signal_filtered = zeros(size(signal));

[b,a] = butter(N,2/Fs*[f_low f_high]);% now let's filter the signal

data_filtered(ind,:) = filtfilt(b,a,data(ind,:));% save as wav file

audiowrite('Immigran_Hiss-filt.wav',data_filtered,Fs);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% display 2B : time / frequency analysis : spectrogram demo

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[sg,fsg,tsg] = specgram(data_filtered(:,1),NFFT,Fs,hanning(NFFT),floor(NFFT*OVERLAP)); % FFT normalisation and conversion amplitude from linear to dB (peak)sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg)); % NB : X=fft(x.*hanning(N))*4/N; % hanning only% apply A weigthing if neededif option_w == 1 pondA_dB = pondA_function(fsg); sg_dBpeak = sg_dBpeak+(pondA_dB*ones(1,size(sg_dBpeak,2))); my_title = ('Spectrogram (dB (A))');else my_title = ('Spectrogram (dB (L))');end% saturation of the dB range : % saturation_dB = 60; % dB range scale (means , the lowest displayed level is XX dB below the max level)min_disp_dB = round(max(max(sg_dBpeak))) - spectrogram_dB_scale;sg_dBpeak(sg_dBpeak<min_disp_dB) = min_disp_dB;% plots spectrogramfigure(3);imagesc(tsg,fsg,sg_dBpeak);colormap('jet');axis('xy');colorbar('vert');gridtitle([my_title ' / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(fsg(2)-fsg(1)) ' Hz ']);xlabel('Time (s)');ylabel('Frequency (Hz)');function pondA_dB = pondA_function(f) % dB (A) weighting curve

n = ((12200^2*f.^4)./((f.^2+20.6^2).*(f.^2+12200^2).*sqrt(f.^2+107.7^2).*sqrt(f.^2+737.9^2))); r = ((12200^2*1000.^4)./((1000.^2+20.6^2).*(1000.^2+12200^2).*sqrt(1000.^2+107.7^2).*sqrt(1000.^2+737.9^2))) * ones(size(f)); pondA = n./r; pondA_dB = 20*log10(pondA(:));endfunction [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)% FFT peak spectrum of signal (example sinus amplitude 1 = 0 dB after fft).

% Linear averaging

% signal - input signal,

% Fs - Sampling frequency (Hz).

% nfft - FFT window size

% Overlap - buffer overlap % (between 0 and 0.95)

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_spectrum = 0; for i=1:spectnum start = (i-1)*offset; sw = signal((1+start):(start+nfft)).*window; fft_spectrum = fft_spectrum + (abs(fft(sw))*4/nfft); % X=fft(x.*hanning(N))*4/N; % hanning only

end fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling

% one sidded fft spectrum % Select first half

if rem(nfft,2) % nfft odd

select = (1:(nfft+1)/2)'; else select = (1:nfft/2+1)'; endfft_spectrum = fft_spectrum(select);freq_vector = (select - 1)*Fs/nfft;end

The passbands and stopbands for the Signal Processing Toolbox filter functions are actually normalised to the closed interval [0, pi] radians. It is calculated as the passband and stopband frequency in Hz divided by the Nyquist frequency.

Example for frequency calculations for the filters your signal:

`Ts = 0.002; % Sampling Interval (seconds)`

Fs = 1/Ts; % Sampling Frequency (Hz)

Fn = Fs/2; % Nyquist Frequency

Note that the highest frequency you can design in your filter is the Nyquist frequency, here 250 Hz.

An example of a Chebyshev Type II filter design for filtering an EKG signal is here:

`Fs = 250; % Sampling Frequency (Hz)Fn = Fs/2; % Nyquist Frequency (Hz)`

Wp = [1.0 100]/Fn; % Passband Frequencies (Normalised)

Ws = [0.5 101]/Fn; % Stopband Frequencies (Normalised)

Rp = 10; % Passband Ripple (dB)

Rs = 50; % Stopband Ripple (dB)

[n,Ws] = cheb2ord(Wp,Ws,Rp,Rs); % Filter Order

[z,p,k] = cheby2(n,Rs,Ws); % Filter Design

[sosbp,gbp] = zp2sos(z,p,k); % Convert To Second-Order-Section For Stability

figure(3)freqz(sosbp, 2^17, Fs)EKG_filt = filtfilt(sosbp, gbp, EKG_original); % Filter EKG Signal

The filter passband goes from 1 Hz to 100 Hz.

## Best Answer