xwvd

Cross Wigner-Ville distribution and cross smoothed pseudo Wigner-Ville distribution

Syntax

d = xwvd(x,y)

d = xwvd(x,y,fs)

d = xwvd(x,y,ts)

d = xwvd(___,"smoothedPseudo")

d = xwvd(___,"smoothedPseudo",twin,fwin)

d = xwvd(___,"smoothedPseudo",NumFrequencyPoints=nf)

d = xwvd(___,MinThreshold=thresh)

[d,f,t] = xwvd(___)

xwvd(___)

Description

d = xwvd(x,y) returns the cross Wigner-Ville distribution of x and y.

d = xwvd(x,y,fs) returns the cross Wigner-Ville distribution when x and y are sampled at a rate fs.

d = xwvd(x,y,ts) returns the cross Wigner-Ville distribution when x and y are sampled with a time interval ts between samples.

d = xwvd(___,"smoothedPseudo") returns the cross smoothed pseudo Wigner-Ville distribution of x and y. The function uses the length of the input signals to choose the lengths of the windows used for time and frequency smoothing. This syntax can include any combination of input arguments from previous syntaxes.

d = xwvd(___,"smoothedPseudo",twin,fwin) specifies the time window twin and the frequency window fwin used for smoothing. To use the default window for either time or frequency smoothing, specify the corresponding argument as empty, [].

d = xwvd(___,"smoothedPseudo",NumFrequencyPoints=nf) computes the cross smoothed pseudo Wigner-Ville distribution using nf frequency points. You can specify twin and fwin in this syntax, or you can omit them.

d = xwvd(___,MinThreshold=thresh) sets to zero those elements of d whose amplitude is less than thresh. This syntax applies to both the cross Wigner-Ville distribution and the cross smoothed pseudo Wigner-Ville distribution.

[d,f,t] = xwvd(___) also returns a vector of frequencies f and a vector of times t at which d is computed.

xwvd(___) with no output arguments plots the real part of the cross Wigner-Ville or cross smoothed pseudo Wigner-Ville distribution in the current figure.

Examples

Cross Wigner-Ville Distribution of Signals

Generate two signals sampled at 1 kHz for 1 second and embedded in white noise. One signal is a sinusoid of frequency 150 Hz. The other signal is a chirp whose frequency varies sinusoidally between 200 Hz and 400 Hz. The noise has a variance of $0.1^2$.

fs = 1000;
t = (0:1/fs:1)';

x = cos(2*pi*t*150) + 0.1*randn(size(t));
y = vco(cos(3*pi*t),[200 400],fs) + 0.1*randn(size(t));

Compute the Wigner-Ville distribution of the sum of the signals.

wvd(x+y,fs)

Compute and plot the cross Wigner-Ville distribution of the signals. The cross-distribution corresponds to the cross-terms of the Wigner-Ville distribution.

xwvd(x,y,fs)

Cross Wigner-Ville Distribution of Chirps

Generate a two-channel signal that consists of two chirps. The signal is sampled at 3 kHz for one second. The first chirp has an initial frequency of 400 Hz and reaches 800 Hz at the end of the sampling. The second chirp starts at 500 Hz and reaches 1000 Hz at the end. The second chirp has twice the amplitude of the first chirp.

fs = 3000;
t = (0:1/fs:1-1/fs)';

x1 = chirp(t,1400,t(end),800);
x2 = 2*chirp(t,200,t(end),1000);

Store the signal as a timetable. Compute and plot the cross Wigner-Ville distribution of the two channels.

xt = timetable(seconds(t),x1,x2);

xwvd(xt(:,1),xt(:,2))

Use Cross Wigner-Ville Distribution to Estimate Instantaneous Frequency

Compute the instantaneous frequency of a signal by using a known reference signal and the cross Wigner-Ville distribution.

Create a reference signal consisting of a Gaussian atom sampled at 1 kHz for 1 second. A Gaussian atom is a sinusoid modulated by a Gaussian. Specify a sinusoid frequency of 50 Hz. The Gaussian is centered at 64 milliseconds and has a variance of $0.{01}^2$.

fs = 1e3;
t = (0:1/fs:1-1/fs)';

mu = 0.064;
sigma = 0.01;
fsin = 50;

xr = exp(-(t-mu).^2/(2*sigma^2)).*sin(2*pi*fsin*t);

Create the "unknown" signal to analyze, consisting of a chirp. The signal starts suddenly at 0.4 second and ends suddenly half a second later. In that lapse, the frequency of the chirp decreases linearly from 400 Hz to 100 Hz.

f0 = 400;
f1 = 100;

xa = zeros(size(t));
xa(t>0.4 & t<=0.9) = chirp((0:1/fs:0.5-1/fs)',f0,0.5,f1);

Create a two-component signal consisting of the sum of the unknown and reference signals. The smoothed pseudo Wigner-Ville distribution of the result provides an "ideal" time-frequency representation.

Compute and display the smoothed pseudo Wigner-Ville distribution.

w = wvd(xa+xr,fs,"smoothedPseudo");

wvd(xa+xr,fs,"smoothedPseudo")

Compute the cross Wigner-Ville distribution of the unknown and reference signals. Take the absolute value of the distribution and set to zero the elements with amplitude less than 10. The cross Wigner-Ville distribution is equal to the cross-terms of the two-component signal.

Plot the real part of the cross Wigner-Ville distribution.

[c,fc,tc] = xwvd(xa,xr,fs);
c = abs(c);
c(c<10) = 0;

xwvd(xa,xr,fs)

Enhance the Wigner-Ville cross-terms by adding the ideal time-frequency representation to the cross Wigner-Ville distribution. The cross-terms of the Wigner-Ville distribution occur halfway between the reference signal and the unknown signal.

d = w + c;

d = abs(real(d));

imagesc(tc,fc,d)
axis xy
colorbar

Identify and plot the high-energy ridge corresponding to the cross-terms. To isolate the ridge, find the time values where the cross-distribution has nonzero energy.

ff = tfridge(c,fc);

tv = sum(c)>0;

ff = ff(tv);
tc = tc(tv);

hold on
plot(tc,ff,"r--",linewidth=2)
hold off

Reconstruct the instantaneous frequency of the unknown signal by using the ridge and the reference function. Plot the instantaneous frequency as a function of time.

tEst = 2*tc - mu;
fEst = 2*ff - fsin;

plot(tEst,fEst)

Input Arguments

x, y — Input signals

vectors | timetables

Input signals, specified as vectors or MATLAB^® timetables each containing a single vector variable. x and y must both be vectors or both be timetables and must have the same length.

• If x and y are timetables, then they must contain increasing finite row times.

• If a timetable has missing or duplicate time points, you can fix it using the tips in Clean Timetable with Missing, Duplicate, or Nonuniform Times.

If the input signals have odd length, the function appends a zero to make the length even.

Example: cos(pi/8*(0:159))'+randn(160,1)/10 specifies a sinusoid embedded in white noise.

Example: timetable(seconds(0:5)',rand(6,1)) specifies a random variable sampled at 1 Hz for 4 seconds.

Data Types: single | double
Complex Number Support: Yes

fs — Sample rate

2*pi (default) | positive numeric scalar

Sample rate, specified as a positive numeric scalar.

ts — Sample time

duration scalar

Sample time, specified as a duration scalar.

twin, fwin — Time and frequency windows

vectors of odd length

Time and frequency windows used for smoothing, specified as vectors of odd length. By default, xwvd uses Kaiser windows with shape factor β = 20.

• The default length of twin is the smallest odd integer greater than or equal to round(length(x)/10).

• The default length of fwin is the smallest odd integer greater than or equal to nf/4.

Each window must have a length smaller than or equal to 2*ceil(length(x)/2).

Example: kaiser(65,0.5) specifies a 65-sample Kaiser window with a shape factor of 0.5.

nf — Number of frequency points

2*ceil(length(x)/2) (default) | integer

Number of frequency points, specified as an integer. This argument controls the degree of oversampling in frequency. The number of frequency points must be at least (length(fwin)+1)/2 and cannot be greater than the default.

thresh — Minimum nonzero value

-Inf (default) | real scalar

Minimum nonzero value, specified as a real scalar. The function sets to zero those elements of d whose amplitudes are less than thresh.

Output Arguments

d — Cross Wigner-Ville distribution

matrix

Cross Wigner-Ville distribution, returned as a matrix. Time increases across the columns of d, and frequency increases down the rows. The matrix is of size N_f × N_t, where N_f is the length of f and N_t is the length of t.

f — Frequencies

vector

Frequencies, returned as a vector.

• If the input has time information, then f contains frequencies expressed in Hz.

• If the input does not have time information, then f contains normalized frequencies expressed in rad/sample.

t — Time instants

vector

Time instants, returned as a vector.

• If the input has time information, then t contains time values expressed in seconds.

• If the input does not have time information, then t contains sample numbers.

The number of time points is fixed as 4*ceil(length(x)/2).

More About

Cross Wigner-Ville Distribution

For continuous signals x(t) and y(t), the cross Wigner-Ville distribution is defined as

$${\mathrm{XWVD}}_{x,y}(t,f)={\displaystyle {\int }_{-\infty }^{\infty }x\left(t+\frac{\tau }{2}\right)y^{*}\left(t-\frac{\tau }{2}\right)e^{-j2\pi f\tau }d\tau .}$$

For a discrete signal with N samples, the distribution becomes

$${\text{XWVD}}_{x,y}(n,k)={\displaystyle \sum _{m=-N}^{N}x(n+m/2)y^{*}(n-m/2)e^{-j2\pi km/N}}.$$

For odd values of m, the definition requires evaluation of the signal at half-integer sample values. It therefore requires interpolation, which makes it necessary to zero-pad the discrete Fourier transform to avoid aliasing.

The cross Wigner-Ville distribution contains interference terms that often complicate its interpretation. To sharpen the distribution, one can filter the definition with lowpass windows. The cross smoothed pseudo Wigner-Ville distribution uses independent windows to smooth in time and frequency:

$${\text{XSPWVD}}_{x,y}^{g,H}(t,f)={\displaystyle {\int }_{-\infty }^{\infty }g(t)H(f)x\left(t+\frac{\tau }{2}\right)y^{*}\left(t-\frac{\tau }{2}\right)e^{-j2\pi f\tau }d\tau .}$$

References

[1] Cohen, Leon. Time-Frequency Analysis: Theory and Applications. Englewood Cliffs, NJ: Prentice-Hall, 1995.

[2] Mallat, Stéphane. A Wavelet Tour of Signal Processing. Second Edition. San Diego, CA: Academic Press, 1999.

[3] Malnar, Damir, Victor Sucic, and Boualem Boashash. "A cross-terms geometry based method for components instantaneous frequency estimation using the cross Wigner-Ville distribution." In 11th International Conference on Information Sciences, Signal Processing and their Applications (ISSPA), pp. 1217–1222. Montréal: IEEE^®, 2012.

Extended Capabilities

C/C++ Code Generation
Generate C and C++ code using MATLAB® Coder™.

Version History

Introduced in R2018b

R2023a: Timetable support for code generation

The xwvd function supports timetable inputs for code generation.