Maximum Likelihood carrier phase estimation of $A\cos(2\pi f_ct)$ This question refers to the example 5.2-1 of the book "Digital Communications, 5th Edition by John Proakis, Masoud Salehi" (example 6.2-1 in the 4th edition of the book) which asks for the estimate carrier phase $\phi$ that maximizes the log-likelihood function and equivalently the likelihood function.

The received signal is:

$r(t)=A\cos(2\pi f_ct+\phi)+n(t)$

where $A$ is the amplitude of the transmitted signal (known), $f_c$ is the carrier frequency (known) and $n(t) $ is the AWGN noise added.

The log-likelihood function is given as:

$\Lambda_L(\phi)=\frac{2A}{N_0}\int_{T_0}r(t)\cos(2\pi f_ct+\phi)dt$

where $\frac{N_0}{2}$ is the power spectral density of AWGN (known constant for all frequencies)
In order to get the maximum value of this we require:

$\frac{d\Lambda_L(\phi)}{d\phi}=0 \\
\Rightarrow \int_{T_0}r(t)\sin(2\pi f_ct+\hat{\phi}_{ML})dt=0 \\
\Rightarrow \hat{\phi}_{ML}=-\tan^{-1}\Bigg[\frac{\int_{T_0}r(t)\sin(2\pi f_ct)dt}{\int_{T_0}r(t)\cos(2\pi f_ct)dt}\Bigg]$

Everything OK until the last line. I don't understand how this equation is solved.

for the last line
$$\sin(2 \pi  f_{c}t + \phi_{ML} ) =\sin(2 \pi  f_{c}t)\cos (\phi_{ML}) +\cos(2 \pi  f_{c}t)\sin(\phi_{ML})
$$
so the integral could be written as:

\begin{align}
 \int_ {T_{0}} r(t)(\sin(2 \pi  f_{c}t)\cos (\phi_{ML}) +\cos(2 \pi  f_{c}t)\sin(\phi_{ML}))&= 0\\
\Rightarrow \cos(\phi_{ML})\int_ {T_{0}} r(t)\sin(2 \pi  f_{c}t) + \sin(\phi_{ML})\int_ {T_{0}} r(t) \cos(2 \pi  f_{c}t) &= 0
\end{align}

Therefore,
$$
\tan(\phi_{ML})= -\frac{\int_ {T_{0}} r(t)\sin(2 \pi  f_{c}t)}{\int_ {T_{0}} r(t) \cos(2 \pi  f_{c}t)}
$$

Hence the result.


How to filter almost periodic noise? I have signal that looks like this:

As you can see there are some almost periodic pulses at the background with amplitude about 1000. I am saying almost, because, if you zoom in, you may see that they slightly differ from each other, plus there may be intervals without them, something like between 35000 and 40000.

How can I get rid of these pulses?

csv with data

I think you have a difficult problem here. Smoothing you signal (for example, with a moving average) probably won't help you. I guess you could try passing your signal through a narrowband notch filter whose notch center frequency is the frequency of your unwanted pulses. MBaz is correct. To think about an appropriate filter means you'd have to know what your signal looks like if it were not contaminated by those unwanted pulses. For example, does your uncontaminated signal look like the following?

I realize I'm not answering your question in any meaningful way, but I was thinking: I assume your noisy signal is the output of an A/D converter sampling an analog voltage signal that's badly contaminated by some sort of electrical noise, some electrical interference. Instead of eliminating the noise in the digital samples I think you should do everything you possibly can to eliminate the electrical noise in your original analog signal.

Dec. 12, 2015 Update based on Vlad's Comment:
@Vladimir: Ah, interesting! What you wrote makes much sense. Those pulses are not noise contamination, they are part of your signal. To attenuate those pulses, three things off the top of my head:

[1] I'd like to see the result of applying your signal to a moving average process.

[2] Although I have no meaningful experience with such a filter, I wonder what would be the result of applying your signal to a median filter.

[3] There's a scheme in the "DSP Tricks" chapter of my DSP book describing how to attenuate very narrow impulsive noise in a signal. I don't know if that scheme would be useful with regard to your pulses but who knows? If you'd like me to send you information on this scheme, send me a private e-mail at: R_dot_Lyons_@_ieee_dot_org.


What features of using Fourier transform method in design FIR filter I want to know when the  Fourier  transform method is best for design  FIR filter

When you know the desired frequency response, and the time shape is irrelevant.


Generalized linear-phase filter I know I can write the frequency response of a Generalized linear-phase system as:
$$H(e^{j\omega}) = A(\omega) e^{-j\left(\alpha \omega - \beta\right)}$$
where $A(\omega)$ is real.

I need to prove that $A(\omega)$ is real and I can't find anywhere the demonstration.

Let's assume that $A(\omega)$ is complex valued:

$$A(\omega)=B(\omega)e^{j\phi_A(\omega)}\tag{1}$$

where $B(\omega)$ is a (possibly bipolar) real-valued function, and $\phi_A(\omega)$ is a (real-valued) phase function.

The phase of

$$H(e^{j\omega})=A(\omega)e^{-j(\alpha\omega-\beta)}=B(\omega)e^{j\phi_A(\omega)}e^{-j(\alpha\omega-\beta)}=B(\omega)e^{j\phi(\omega)}\tag{2}$$

is then

$$\phi(\omega)=\phi_A(\omega)-(\alpha\omega-\beta)\tag{3}$$

(where I define the phase such that jumps by $\pi$ are taken care of by the sign of $B(\omega)$, not by $\phi(\omega)$). For the phase $\phi(\omega)$ to be linear, $\phi_A(\omega)$ must be constant or linear:

$$\phi_A(\omega)=a\omega+b\tag{4}$$

With ($4$) the total phase becomes

$$\phi(\omega)=(a-\alpha)\omega+(b-\beta)=-\alpha'\omega+\beta'\tag{5}$$

Combining $(2)$ and $(5)$ we get

$$H(e^{j\omega})=B(\omega)e^{j\phi(\omega)}=B(\omega)e^{-j(\alpha'\omega-\beta')}\tag{6}$$

which has the same form as the equation in your question. Note that in $(1)$ we have defined $B(\omega)$ to be real-valued.

This shows that $(6)$ (or, equivalently, the equation in your question) is the most general form of the frequency response of a linear phase filter. If $A(\omega)$ were chosen to be complex-valued, then $A(\omega)$ must have a linear phase (otherwise $H(e^{j\omega})$ wouldn't have a linear phase), and that linear phase can be absorbed by the constants $\alpha$ and $\beta$, such that the expression for the frequency response can always be reduced to the form $(6)$ with a real-valued function $B(\omega)$.


Fourier Transform of the full Morlet wavelet In 2014 someone asked here the Fourier transform of the Morlet wavelet; link below:

Fourier Transform of Morlet wavelet Function?

However, it was the approximated Morlet wavelet not written with the canonical Gaussian function.

Can somebody help to find the Fourier Transform of the full scaled Morlet wavelet below:

$$\psi_{\mu}(\tfrac{t}{\mu}) = \frac{1}{\mu} \times \frac{1}{\sqrt{2\pi} \, \sigma} \, e^{\frac{ -\left(\frac{t}{\mu} \right)^2}{2\sigma^2}} \big[e^{j\omega_c \frac{t}{\mu}}-e^{-\frac{1}{2} \omega_c^2}\big]$$

where $\omega_c = 2\pi f_c$ and $f_c$ is the centre or carrier frequency of the wavelet (the frequency to which the wavelet oscillates in the temporal domain). Also, $\mu$ is the scaling factor defined as $\mu = \frac{\omega_c}{\omega_{\mu}}$, where $\omega_{\mu}$ is the analysing frequency. Indeed, when $\omega_{\mu} = \omega_{c} \Rightarrow \mu = 1$ we have the mother Morlet wavelet.

In addition, the second term in the brackets is the correction term necessary to enforce zero mean to the Morlet wavelet.

I'd much appreciate a detailed step-by-step demonstration if possible to check where I'm getting stuck in my own demonstration. It could be done in a paper if quicker and a photo sent to me; I'll then publish here the full demonstration here.

Thanks in advance for any help.

if you wanna do this with ordinary frequency $f$ (rather than angular frequency $\omega$) youcan start with

$$ \mathscr{F}\Big\{ x(t) \Big\} \triangleq \int\limits_{-\infty}^{\infty} x(t) \, e^{-j 2 \pi f t} \, dt $$

$$ \mathscr{F}\left\{ e^{-\pi t^2} \right\} = e^{-\pi f^2} $$

$$ \mathscr{F}\left\{ e^{-\pi \alpha t^2} \right\} = \frac{1}{\sqrt{\alpha}} e^{-\frac{\pi}{\alpha} f^2} \qquad \alpha > 0 $$

$$ \mathscr{F}\left\{ e^{-\pi \alpha (t-\tau)^2} \right\} = \frac{1}{\sqrt{\alpha}} e^{-\frac{\pi}{\alpha} f^2} e^{-j 2 \pi f \tau} $$


Why does it make sense to have noise of a different shape than the state vector? Usually when I see the Kalman filter in lectures (or other resources like 1), the state equation is

$$\mathbf{x}_{k+1} = A_k \mathbf{x}_k + B_k a_k + r_k^{(s)},$$

and the measurement equation is

$$z_k = H \cdot \mathbf{x}_k + r_k^{(m)}$$

where $r_k^{(m)}, r_k^{(s)}$ is normal distributed noise of the same shape as the state vector $\mathbb{x}_k$.

Then, the covariance in the prediction step is usually stated as

$$P_{k+1}^{(P)} = A P_k A^T + C_k^{(r_s)}$$

where $C_k^{(r_s)}$ is the covariance matrix of the system noise.

However, in one lecture the state equation was stated as

$$\mathbf{x}_{k+1} = A_k \mathbf{x}_k + B_k a_k + G r_k^{(s)}$$
and the prediction for the covariance matrix of the system noise was
$$P_{k+1}^{(P)} = A P_k A^T + G Q G^T$$

Why does it make sense to model the noise as something which gets multiplied by a matrix $G$ in contrast to simply having a noise vector which has the same shape as the state vector?

(See my article for detailed explanations of the variables / how I understand the Kalman filter right now.)

Why does it make sense to model the noise as something which gets multiplied by a matrix $G$ in contrast to simply having a noise vector which has the same shape as the state vector?

The usual reason is because the model has some states that are not directly driven by the process noise. For example, suppose your states are:
$$
\mathbf{x}_k = \left [ \begin{array}{c}
d_k\\
v_k\\
a_k
\end{array} \right]
$$
where $d_k$ is displacement, $v_k$ is velocity, and $a_k$ is acceleration all at time index $k$.

Suppose, as with a car, you can only change the acceleration component of the state (by pressing on the accelerator or the brake [decelerator]):

$$
\mathbf{x}_{k+1} = \mathbf{A} \mathbf{x}_{k} +
\left [ \begin{array}{c}
0\\
0\\
1
\end{array} \right] v_k
$$

In this case,
$$
G = \left [ \begin{array}{c}
0\\
0\\
1
\end{array} \right]
$$


best technique to remove noise from bio signal data I'm still learning about signal processing so sorry if this is an intro problem. I want to remove noise from the following data so I can separate each wave and run analysis on it.

The data is sampled at 41 hz, the red is the raw data and the blue is the data that I applied a lowess filter to with a window size of 41. As you can see it is smoothed to some extent but there is still a lot of peaks and noise in the mix.

I read that you can use RNN or PCA or FFTs to remove noise for bio-signal data, in this case, what are the best techniques to remove noise from the red data so that the blue line is just 5 waves that is easily identifiable?

My ideal result is that after applying the filter I can easily tell that there is 5 waves in this period of data and the resulting index of said waves.
Any input is appreciated!

As a first appraoch I would suggest to use a mean filter, which will calculate the average of the last n values.

$$y[n] = 1/N \sum_{k=0}^{N-1} x[n-k]$$

If you have access to the full signal (offline processing) you can perform forward and backward filtering which will avoid the delay introduced by a single filter. (Check matlab's filtfilt function for example)


Concept of finding forced response from difference equation? I know how the process works, first from difference equation make a homogeneous equation and then find natural response , then find a particular response (in the form of input) substitute it back in the equation and find constants.  And we get a forced response.  I have a doubt in finding constants.
This example is Drill problem 2.12 in "Signals and Systems" by Simon Haykin:

For the equation:
$$y[n] - (0.25)*y[n-2] = 2*x[n]+x[n-1]$$
find a forced response for $x[n]=u[n].$

I have found the natural response and for particular response. Say, $y[n] = k*u[n]$
Now the equation is:
$$k*u[n] - 0.25*k*u[n-2] = 2*u[n] + u[n-1]$$
The answer according to the book is $k=4$, but that is true only for $n>=2$. Can somebody explain what happens if $2>n>0$; $k=4$ does not seem correct in this case, as $u[n-2]=0$ for $n<2$ and $u[n-1]=0$ for $n<1$.
Why can we just substitute all $u[n-c] = 1$, in the difference equation for finding value of constants in particular solution?

considering your equation,

$$k\times u[n] - 0.25\times k\times u[n-2] = 2\times u[n] + u[n-1]$$

the left half of the equation will result in a series as shown below:

$$\begin{array}{cccccccc} [k & k & \frac{3}{4}k & \frac{3}{4}k & \frac{3}{4}k & \frac{3}{4}k  & \ldots &]\end{array}$$

the right hand of the equation is
$$\begin{array}{cccccccc} [2 & 3 & 3 & 3 & 3 & 3 & \ldots & ]\end{array}$$

for these two series to be equal,
$$\begin{array}{cccccccc} k = 2 \\ k=3\\ \frac{3}{4}k=3 \\ \frac{3}{4}k=3 \\ \vdots & \end{array}$$

It is not possible to equate these two sequences without removing first two points. In that case, $k = 4$. The sequences are equal or the equation holds for $k = 4$ and $n \ge 2$.

What we are getting here is total solution, so it will have the natural response and forced response. The output get stabilized and reach the stable value only after $n \ge 2$. This 'transients' before $n < 2$ are due to the natural response.


Derivation of Toeplitz Matrix I'm having a difficult time understanding why the matrix for LTI systems is a Toeplitz matrix. I can see why $h_{n,m} = h_{n' + q,m' + q}$ given that $n' = n - q$ and $m' = m - q$, and $$\sum_{m'= -\infty}^\infty h_{n' + q,m' + q} x[m'] = y[n']$$ and $$\sum_{m= -\infty}^\infty h_{n,m} x[m] = y[n],$$ but the leap to the equivalence $h_{n,m} = h_{n + q,m + q}$ escapes me.

This reference should help.  If you carry out the matrix-vector multiplication, you will see they are the same.

Also note that it is only possible to express the convolution as a matrix multiplication when the filter impulse response is of finite duration, i.e. a FIR filter.

The equivalence $h_{n,m} = h_{n+q,m+q}$ comes from the "TI" part of LTI, that is, time-invariance. If you delay the input by $q$ samples, the output should be the same, just delayed in time by $q$.


Which type of fliter should be designed to remove 50 Hz Actually I am using DSP kit TMS320C6711 through which I want to remove 50Hz component out of the signal. I am unable to interface Code Composer studio (v2) with matlab. Therefore, I have designed FIR notch filter using Code composer Studio in which coefficients are taken from fdatool of matlab. Kaiser window has been used for filtering. The problem I am facing here is when lower and upper cutoff frequencies are 35Hz and 65Hz respectively only I am able to manage significant attenuation. When I am trying to get the cutoff frequencies closer, attenuation level is reducing to unsatisfactory levels.
So, please help me to decide the window to be used and order to be selected so that I could achieve significant attenuation of 50Hz component with sharp cutoff frequencies.

Highly selective notch filters require a rather long impulse response, which makes FIR an inefficient approach. A simple IIR notch filter is easily designed by zero/pole placement:

Take a 2nd order recursive filter and place two zeros at $\exp(i \omega_0)$ and $\exp(-i \omega_0)$ on the unit circle, where $\omega_0=2\pi\frac{f}{f_s}$ is the normalised center frequency of your notch. Then place two poles just inside the unit circle behind the zeros, namely at $r\exp(i\omega_0)$ and $r\exp(-i\omega_0)$ with $0<r<1$. These two poles will cancel the effect of the zeros beyond a certain distance. The closer the poles are to the zeros (i.e. the closer $r$ is to $1$), the sharper the notch will be. Picking the right $r$ therefore controls your bandwidth but also gives you the ability to control numerical stability and error accumulation for limited precision implementations. If you more the poles too close to the unit circle, the filter will start to behave weirdly.

Using this recipe to calculate your filter coefficients we can directly start with the z-domain transfer function
$$H(z) = \frac{\left(z-\exp(i\omega_0)\right)\left(z-\exp(-i\omega_0)\right)}{\left(z-r\exp(i\omega_0)\right)\left(z-r\exp(-i\omega_0)\right)}$$
which can be expanded to
$$H(z)=\frac{z^2-z\left(\exp(i\omega_0)+\exp(-i\omega_0)\right)+1}{z^2-z r\left(\exp(i\omega_0)+\exp(-i\omega_0) \right)+r^2}=\frac{z^2-2\cos(\omega_0)z+1}{z^2-2r\cos(\omega_0)z+r^2}$$
and cancelling $z^2$ gives you the transfer function in terms of unit delays:
$$H(z) = \frac{1-2\cos(\omega_0)z^{-1}+z^{-2}}{1-2r\cos(\omega_0)z^{-1}+r^2z^{-2}} $$
Using the definition of the transfer function $H(z)=\frac{Y(z)}{X(z)}$ for input and output signal z-transforms $X(z)$ and $Y(z)$ respectively results in the z-transform of the filter recursion:
$$Y(z)\left(1-2r\cos(\omega_0)z^{-1}+r^2z^{-2} \right)=X(z)\left(1-2\cos(\omega_0)z^{-1}+z^{-2} \right)$$
Inverting the z-transform takes us to discrete time domain:
$$y[n]-2r\cos(\omega_0)y[n-1]+r^2y[n-2]=x[n]-2\cos(\omega_0)x[n-1]+x[n-2]$$
Finally, we can solve for the output signal at time $n$ and get the recursive filter equation
$$y[n]=x[n]-2\cos(\omega_0)x[n-1]+x[n-2]+2r\cos(\omega_0)y[n-1]-r^2y[n-2]$$
This equation gives you the most recent output sample in terms of the three most recent input samples and the prior two output samples. As initial condition you can assume that all recent samples are 0.


When to use Rician fading vs Gain for Directional Antenna This is a pretty simple question. [The answer might not be simple :( ]

If I am using a directional antenna, can I used Antenna gain $G_T=1$ as I am considering Rician Fading between transmitter and receiver?

On the contrary, if I am using $G_T>1$, then do I use Rayleigh fading instead of Rician?

Explanations of System Model:

The antenna is mounted under the roof-top in an indoor environment. (e.g. indoor stadiums, warehouses, office building/rooms etc.). The Antenna is directed towards the ground, so that the signal-loss in unwanted direction is minimized, (unlike isotropic antenna where energy is radiated in $360^o$).

The gain is given by the antenna used, and for directional antennas the gain would be much greater than 1. $G_T$ is typically given as the antenna gain compared to an isotropic antenna which radiates or receives equally in all directions.  (So $G_T=1$ corresponds to a perfect isotropic antenna, not a directional one).

For example, the azimuth pattern below is for a Yagi antenna with a gain of 15 dBi, which shows the highest gain is achieved when the antenna is pointed in the 90° direction, along with other sidelobes that are also sensitive but significantly less gain (for example approximately -4 dBi at 45° and 135°). Not shown but the elevation pattern also has sidelobes that approach 0 dBi.

image source: Cisco white-paper on antenna patters

Rayleigh fading is a statistical model for the channel with multiple reflected propagation paths without any particular path being overly dominant, while Ricean results is the model for the channel that has multiple reflected propagation paths along with a dominant (usually direct) path. Rayleigh fading is typically a reasonable model for heavily urban areas while Rician is often used for SatCom. That said, a highly directional antenna can have the benefit of minimizing reflected paths from other spatial directions and therefore can make an otherwise Rayleigh channel appear to be more Rician, assuming that the antenna is pointed at a propagation path that is significantly strong enough, and the environment is not one that would still result in multiple paths of random delay within the main lobe of the antenna without a single dominant propagation path.

So a highly directional antenna can make a Rayleigh channel into a Rician channel, but that is not guaranteed to do so. It may seem strange at first to state that the antenna can change the channel, but this is really no different than all the buildings, mountains and other objects causing the reflections that also create the channel: the presence or absence of reflections is what creates the channel and the type of antenna is just one of the objects considered in determining the reflection environment.

For further discussion on Rayleigh vs Rician fading see Rayleigh fading with frequency selective fading channel


Steganography - What techniques exist to hide an image in an audio file? Can anyone mention me techniques to hide an image in an audio file ?

Thanks

i dunno if this is a technique oft discussed in the field of steganography, but this is the first technique mentioned for digital data:

you embed your hidden message in the least-significant bit (LSB) of every sample word in the audio.  i have used this to hide (from audibility) metadata about the audio.  i have, in the past, embedded pitch detection data in the LSB of a stream of 16-bit audio words.  the result is that the true audio was really 15-bit and the hidden embeded data (from the POV of audio was random or noisy bits in the LSB) would just sound like a form of added noise, similar to quantization noise (as if the LSB was always zero).  to best mask this metadata from the sound, noise-shaping or error-shaping is employed at the 15 bit quantization level.

the result is that your meta-data is a stream of bits, one bit per sample of audio.  you'll need to worry about framing and synchronization of that data so you know where the beginning of every meta-data frame is.  the way i did it was similar to the old SDLC in that i chose some "magic number" as a prefix (this is in the serial data from the LSBs of the audio words) to mark the beginning of a frame of metadata.  say that prefix is 0b11111111, eight 1 bits in a row.  if the meta data ever needed to have seven 1 bits in a row, whether or not the next bit was a 1, i would insert an extra 0 into the stream that the receiver knew to take out.  each time this would occur in the serial data, this would make the serial stream 1 bit longer for the frame than would otherwise be required.  that way eight 1 bits in a row could only mean a frame start marker.

actually i would not use 0b11111111 because it would occur in data often, so i analyzed a stream of actual metadata with a C program to see what 8-bit sequences occurred least often.  and of the least common 8-bit sequences i would pick the one with the first 7 bits occurring least often.  that would be my choice for a frame marker.  whatever was the last bit of that 8-bit marker, i would insert the opposite of that bit if i ever saw the uncommon 7-bit sequence in the stream of data and when the receiver would see the same 7 bits (and it wasn't the 8-bit marker), it would remove that bit from the stream.  so there was no possibility of an accidental combination of data falsely screwing up my metadata framing.

so whether this metadata is an image or something else doesn't matter.  it's just data.  of course, if the audio is compressed or processed in any normal manner, it will kill this metadata unless you strip it out and put it back in after processing.


Mathematically speaking, is a "signal" a function or the set of outputs from a function In engineering it is well known that we abuse the notation of a function $f: A \to B$ with the image of that function $f(t) = \left\{y \in B\mid f(t) = y, t \in A\right\}$

When we say: let $f(t) = \sin(\omega t)$ be a signal, do we mean a function or the set of output from that function i.e. is signal a function or the output?
Similarly, looking at the diagram below, would one say that $\epsilon(t)$ is a set of inputs (i.e. scalars) going into the controller, or a function?

A signal is a physical quantity (e.g. voltage) carrying information, or a set of values (e.g. samples in discrete case) of the given function for different values of the underlying independent variable.

In the diagram above $\epsilon(t)$ is a mathematical modeling (i.e. function) that characterizes the set of inputs (i.e. scalers) going into the controller.

Strictly speaking, a signal is the "observed" (measurable) physical quantity while a function is a mathematical model (abstract representation) describing the relationship between the involved variables.


Spectrograms for neural nets General Question

Given an audio file, say a 16-bit wav, what are some standard methods to preprocess a spectrogram of this wav so that it may be fed into a neural net?

Context

In Lee et al's 2009 NIPS paper on applying deep belief networks (DBNs) to speech recognition, they briefly mention in Sections 2.2 and 3.1 preprocessing a spectrogram using PCA whitening.

NB: I don't have image-processing background as applied to machine learning...

Question

How may we reproduce the training data that is fed into the DBN in the above NIPS article, section 4.1 for speaker identification in the TIMIT dataset?

"We trained the CDBN on unlabeled TIMIT data (as described in Section 3.1)...We computed the first and second-layer CDBN features using the
spectrogram as input."

Section 3.1 mentinons 20ms window sizes with 10ms overlaps to generate the spectrograms and then applying PCA whitenting. Does this spectrogram -> whitening process produce an $n \times  m$ or $n \times m \times t$ array for each speaker? Or do we get an "image" that is somehow fed directly into a (C)DBN?

Outlook

Based on the NIPS paper, such preprocessing seems standard but I'm (obviously) new to using spectrograms in neural networks.

Any python references on spectrograms are appreciated as related to these questions, for example with librosa, pylab, or this tutorial.

After generating a spectrogram, perhaps the general steps in this PCA whitening overview can be applied to produce relevant input to a DBN?

Example

For many supervised classification tasks, we may fit a model via fit(X, y) where $X$ is an $n \times k$ matrix of $n$ samples with $k$ features and $y$ is an $n \times 1$ vector where each row is a corresponding label.

Using librosa's melspectogram method on an audio file, we produce an $m \times t$ matrix. With, say, $n$ labeled audio files we would have $n$ training samples each of shape $m \times t$ resulting in an $n \times m \times t$ matrix $X$ (and $n \times 1$ matrix $y$). How may we feed this $X$ into a DBN or neural network in general?

The "dimensions" of the spectrogram are not chosen based on where will the spectrogram be fed to but rather depend on your application. Therefore, it is key to understand the spectrogram itself first, as a means of generating features for one or more signals and to an extent, understand the Discrete Fourier Transform (DFT) as well, which is the key operation the spectrogram is based on. In fact, the "dimensions" you are referring to stem entirely from the way the DFT works and determine its spectral "resolution".

If you were to run DFT on some discrete time signal ($x[n]$), of maximum length $N$, you would get back a set of coefficients $X[k]$ that decompose the signal to a weighted sum of sinusoids at different frequencies.

The range of physical frequencies in $x[n]$ is $\left[ 0 \ldots \frac{F_s}{2} \right]$. However, the DFT decomposes the signal into $N_{FFT}$ frequency bins (indexed by $k$) with $N_{FFT} \le N  $. Therefore, the physical frequency that is associated with the $k^{th}$ frequency bin is $f_k = F_s \cdot \frac{k}{N_{FFT}}$.

So what?! So, if you are creating a classifier that operates on some $x_m[n]$ sample obtained from the $m^{th}$ speaker and is trying to figure out the gender of the speaker via the range of frequencies audible, then you would have to tune your $N_{FFT}$ in such a way as to be able to assign different frequencies between 85Hz to 225Hz to different frequency bins. If you did not do that, then your spectral features (the set of $k$ frequency bins) would be insensitive to your class $C$ (here $C \in \left\{ male,female \right\}$ and the classification would report poor performance.

Alright, so, we keep this consideration of spectral resolution in mind. This actually considers the vertical dimension of the spectrogram. Let's have a look at the spectrogram itself now.

The spectrogram is the repetitive application of Short Time Discrete Fourier Transforms (STFT) to a signal in overlapping sections. Let's denote the overlap here with some factor $q \in \left\{ 0 \ldots 1\right)$. When $q=0$ there is no overlap and the spectrogram is evaluated independently on a set of windows (for example $n_{w1} \in \left\{0 \ldots 10 \right\}, n_{w2} \in \left\{11 \ldots 20 \right\}$ and so on). As $q$ approaches $1$, the STFT starts running on overlapping windows (for example $n_{w1} \in \left\{0 \ldots 10 \right\}, n_{w2} \in \left\{5 \ldots 25 \right\}$ and so on).

This overlap is desired in the STFT for a number of reasons that are related to the maths of the FT but also the signal itself and the desired analysis on it. The bottom line is that now we have to worry about one more parameter, $q$.

So what?!? So, to extend the previous example: Before, we were trying to understand what is the gender of the speaker just by examining a sample of their voice. Let's now assume that we have the recording of a conversation between males and females and we want to produce a list of time instances and their gender classification to examine (for example) if males or females tend to dominate a given discussion. Obviously, we have to tune $N_{FFT}$ from above but now we also have to tune $q$. Too much overlap means a high temporal resolution but also large amount of data (and processing). Too little overlap and we might miss some transition as some speaker was interrupted by some other speaker and then they resumed. If this phenomenon happens within the duration of a window, the STFT cannot resolve different males or females. It's a gray area that contains a mixture of both fundamental frequencies and it will produce some decision that might not be accurate.

As you can see, our discussion so far has not involved the classifier itself. The classifier is indeed some $C = f(x_m,\Theta)$ that takes the $m^{th}$ sample and maps it to one of the elements of the set of classes $C$ based on some parameter set $\Theta$. But the only thing we have been talking about so far is how do you construct this set of features that are derived from $x_m$ so that then the $f$ can do its mapping.

Let's take a look at that now: Remember how playing around with $q$ produced overlapping windows? Let's take this to the extreme. Let's think about what happens when $q=0.95$. Then, two neighbouring windows actually look very much alike. They differ only in 5% of the samples! Consequently, their spectra (adjacent vertical lines of the spectrogram) are very much alike. Therefore, it is as if we are feeding the classifier redundant data. We are showing the classifier lots of similar examples. Therefore, continuing with the above example, the classifier starts forming an opinion that MALE means all of these similar data that we are feeding into it when MALE, in our world, in the physical world, is conceptually all the different MALE voices across all recordings. A large STFT overlap could lead to overfitting!

To brake this similarity of adjacent samples, we introduce whitening. What does whitening do? It shuffles the input matrix, to reduce the covariance between adjacent samples, this is all that the whitening matrix $W$ does (but does it in a mathematical way that takes into account the covariance of the input. It produces an optimal shuffling, so optimal that any covariance between samples has been reduced to zero.)

Of course, when you bring Principal Component Analysis (PCA) into the mix, the game is changing slightly because PCA is still on the feature generation side. Applied to our previous example, PCA would enable us to get an indication of how many distinct speakers there might be in a given conversation and express the rest of the samples as coefficients to these distinct speakers (or waveforms).

So:

  As a multi-dimensional array what dimensions should an audio sample's spectrogram have before PCA whitening;

It depends on your application and specifically the spectral ($N_{FFT}$) and temporal ($q$) resolutions necessary to resolve unique examples for the classifier.

  What dimensions should the pre-preocessed object have so that it may be fed into a DBN?

It doesn't matter. You can choose to feed the whole thing and let the network "choose" the best input vectors by adjusting its weights OR you could run feature selection and reduce the dimensionality of your input (and the effort to train the network and the time it takes to train it, etc).

Hope this helps.

EDIT:
The paper states "For the application of CDBNs to audio data, we first convert time-domain signals into spectrograms.
However, the dimensionality of the spectrograms is large (e.g., 160 channels). We apply PCA whitening to the spectrograms and create lower dimensional representations. Thus, the data we feed into the CDBN consists of $n_c$ channels of one-dimensional vectors of length $n_V$ , where $n_c$ is
the number of PCA components in our representation.". Therefore, the input can be considered as "2D".


Interpreting the inverse fourier transform from a graph I'm given a graph of the fourier transform of some function $x(t)$.  The graph is labelled $F(X(\frac{\omega}{\pi}))$ on the y-axis and $\frac{\omega}{\pi}$ on the x-axis.  The graph is plotted only on the interval from 0 to 1 (corresponding to $\omega$ values ranging from 0 to pi).  There are points at

$\left(\frac{\omega}{\pi},F(\frac{\omega}{\pi})\right)  = (0.2,1); (0.4,2); (0.6,3); (0.8,3);$ and  $(1,2)$

and is zero everywhere else.

I'm confident that these nonzero points correspond to a linear combination of unit sample delta functions centered at $\omega= \pi/5, 2\pi/5, 3\pi/5, 4\pi/5$, and $\pi$, respectively.  Since the delta functions are purely real valued, I conjecture that the inverse fourier transform of this function should be a linear combination of cosine functions.  Based upon a table of fourier transforms, I surmise that $$x(t)=\frac{1}{\pi}cos\left( \frac{\pi t}{5}\right) + \frac{2}{\pi}cos\left( \frac{2\pi t}{5}\right) + \frac{3}{\pi}cos\left( \frac{3\pi t}{5}\right) + \frac{3}{\pi}cos\left( \frac{4\pi t}{5}\right) + \frac{2}{\pi}cos\left(\pi t\right)$$.

Have I made the appropriate conclusions?

If there are no negative frequencies, $x(t)$ would have an imaginary component (see http://en.wikipedia.org/wiki/Analytic_signal), but we can just assume that the graph was truncated for convenience. So I think you have got it right.

Note: previously I wrote in this answer that the coefficients should have been multiples $\frac{1}{2\pi}$ but they are correct if we assume the negative frequencies mirror the positive.


Given a system with Transfer Function and its desired output. Is it possible to find the required Input? I'm currently working on an DC servomotor. I managed to find it transfer function using System Identification theory and Matlab.

The transfer function is given as
$$H(Z) = \frac{1 - 0.4952z^{-1} + 0.01826z^{-2}}{1.612z^{-1} - 0.02492z^{-2}}$$

Now, I have a pre-defined velocity profile. That mean at each sampling time $nT_s$, I have to reached a certain speed. In other words, I have a matrix of required output at time $nT_s$.

Now my question is, is it possible to calculate the input to the transfer function $u(t)$ so that the output approximates the desired velocity? If it is possible, could you kindly give me some broad instruction or some reference so that I can study on it? If no, is there any possible theory that can solve my problem, that is find required Input to get desired Output?

My process so far: Well, I simply inverse that transfer function and ask Matlab to find the required input for me. However, the inversed is not proper and thus it is not possible to do that. So I add some fast-orders (first time hearing this, not sure what it means) to the denominator to make it proper. The function runs, but the results is quite... not right.

This is my desired Output vRequired, stored in a $1$x$714$ matrix and plot for clarity.

And this is the result with the function lsim(reservedTf, vRequired):

where reservedTf is my transfer function reserved with fast-orders:

$$  H^{-1}(z) = \frac{z^{-1} - 0.4952 z^{-2} + 0.01826 z^{-3}}{1.612*10^{-10} + 3.224*10^{-5} z^{-1} + 1.612 z^{-2} - 0.02492 z^{-3}} $$

I very much appreciate all your help. Please tell me if there is anything wrong with my post here.

The transfer function is given as ...

Your transfer functions seems wrong. It's not "proper", i.e. there ought to be a $1$ in the denominator. You can try to rewrite it this way, but it would make it non-causal. Not sure if that's a typo or an error of your System Identification

  However, the inversed is not proper

No, your original transfer function is improper, the inverse is just fine.

  If it is possible, could you kindly give me some broad instruction

Make sure your original transfer function is proper, fix this first
Check if your corrected transfer function is causal, stable & minimum phase, i.e. all poles and zeros inside the unit circle. If that's true, you can simply invert it by exchanging poles and zeros.
If you have zeros outside of the unit circle, you can still invert this, but the result will be non-causal. You can still approximate an inverse by either adding bulk delay (if your application can tolerate that) or by least square error time domain solution.


FSK modulation in baseband In a $4$-FSK we have
$$s_i(t)=\cos(2 \pi f_i t)=\cos(2 \pi (f_c+\alpha_i \Delta f)t)$$
with $\alpha_i=2i-5$, and $\Delta_f=\frac 1T$, where $T$ is the time symbol and $i \in [1,4]$.
But with the expression of $s_i$, I have a bandpass signal. To get the baseband signal: $g_i(t)=e^{j2 \pi f_i t}$.
Is it correct?

The expression $g_i(t) = e^{j2\pi f_i t}$ is not the baseband signal--- it is simply the positive frequency portion of your bandpass signal. To convert this to a baseband signal we need to take the further step of frequency translating the modulated signal centered about the carrier frequency $f_c$ to baseband $f=0$ as follows:
$$f_b = g_i(t)e^{-j 2\pi f_c t} = e^{j2\pi f_i t}e^{- j 2\pi f_c t} = e^{j2\pi (f_i-f_c)t}$$
$$= e^{j 2 \pi (f_c+\alpha_i \Delta f - f_c) t)} = e^{j2\pi \alpha_i \Delta ft}$$
This is further demonstrated with the spectral plots below showing the difference between a passband and baseband signal:


Lane detection and tracking using B-snake? I am truly sorry if my question sounds little bit stupid, but I am getting stuck at this part for a long time. I am currently working for my final year project of lane detection for curved road using B-snake. I am using Matlab software, and I have done the image pre processing part by converting the RGB value to grayscale. After the converting I simply use the Canny edge detection but it seems that there is too many noise.

This is the first part that I obtained :

Then, I tried the region props function for which the Matlab code I obtained from this website. The noise has reduced a little bit.

The only region that I want is the road region. But I don't manage to get the ROI that I need. So I just proceed to Hough Transform to obtain the straight line. I try to run the coding of Hough Transform that I found on the internet, but the result is so confusing. I don't know how to change the code so that my result will turn out better.

The image after Hough Transform is applied :

This is my full code:

I=imread('C:\Users\LENOVO\Documents\MATLAB\PSM2\road\draw\test1.jpg'); %         Input image file here
dim = size(I); % Get size of the image and store in an array called 'dim'.%size array matrix; 1 is height, 2 is width

I2 = rgb2gray(I); % Convert to grayscale
BW = edge(I2,'canny',[]); % Edge detection using Canny
cc = bwconncomp(BW);
I4 = labelmatrix(cc);

a_rp = regionprops(cc,'Area','MajorAxisLength','MinorAxislength','Orientation','PixelList','Eccentricity');
idx = ([a_rp.Eccentricity] > 0.99 & [a_rp.Area] > 100 & [a_rp.Orientation] < 60 & [a_rp.Orientation] > -90);

BW2 = ismember(I4,find(idx));
[H,T,R] = hough(BW2);

%  figure, imshow(H,[], 'XData', T, 'YData', R, 'InitialMagnification', 'fit');
%  xlabel('\theta'), ylabel('\rho');
axis on, axis normal, hold on;
P  = houghpeaks(H,50,'threshold',ceil(0.1*max(H(:))));

% Set houghpeaks parameters, threshold unsure
x = T(P(:,2));
y = R(P(:,1));
plot(x,y,'s','color','white');

% Apply median filtering
I3 = medfilt2(I2);

% Find lines and plot them
lines = houghlines(BW,T,R,P,'FillGap',20,'MinLength',10);
figure, imshow(I3),imagesc(I3), hold on
max_len = 0;

for k = 1:length(lines)
    xy = [lines(k).point1; lines(k).point2];
    plot(xy(:,1),xy(:,2),'LineWidth',2,'Color','green');

    % plot beginnings and ends of lines
    plot(xy(1,1),xy(1,2),'x','LineWidth',2,'Color','yellow');
    plot(xy(2,1),xy(2,2),'x','LineWidth',2,'Color','red');
end

showlines = struct(lines);
cellData = struct2cell(showlines);

% X-coordinates are for width
% Y-coordinates are for height
%point1(x y) etc

for i = 1:280
    % 'A' stores all 'x' coordinates of point 1
    A([i,i+1])= [cellData{1,i}];
    % 'B' stores all 'x' coordinates of point 2
    B([i,i+1])= [cellData{2,i}];
    % 'C' stores all 'y' coordinates of point 1
    C([i,i])= [cellData{1,i}];
    % 'D' stores all 'y' coordinates of point 2
    D([i,i])= [cellData{2,i}];
end

I do really hope someone can help me. At least general idea for me to proceed to the next step is enough. Thank you in advance.

Following my comment I suggest you a "dumb" heuristic you could try : Consider your images row-wise. Extract a few features from each row, these could be average RGB value, RGB covariance matrix, entropy etc ... Use a simple clustering algorithm to separate your image into 2 parts : the upper one and the lower one. Let's assume the heuristic has been successful and the lower part is mostly composed of road. We'll just consider this part and run a clustering algorithm on the RGB values. We keep the number of cluster low and we want to take the cluster which is the closest from grey as the road cluster. Once you have a not too bad road segmentation, play around a bit with edge detection to extract only the white marking.


Is Bluetooth Low Energy jamming possible with an SDR like the HackRF on GNURadio? It seems like BLE jamming is technically possible by just jamming the 3 advertising channels 37, 38 and 39 https://github.com/lws803/BLE-jammer and it has been proven to work on an STM32 + 3x NRF24l01 transceivers. Can the hackRF switch between these 3 frequencies fast enough to produce a consistent jamming effect?
Disclaimer: I understand that this is illegal and I would like to know about this for educational purposes only.

I am not revealing any big secrets here on jamming and anti-jamming techniques, nor would I condone creating any such interference. What I am about to say is quite simplistic and well known, but knowing more details in how jamming can take place and being more educated on it in general can help good actors in minimizing vulnerabilities in future designs.
Yes this is completely feasible but I would do this much more simply:
Create a single sinusoidal tone at 15 MHz and at 39 MHz.
Upconvert this tone using a carrier at 2441 MHz.
This will produce two symmetric sidebands with a tone at:
2402 MHz
2426 MHz
2456 MHz
2480 MHz
So there is one additional channel rather than just the targeted three channels in the hopping approach, but still significantly better and more effective over trying to jam the entire spectrum from 2402 to 2420 MHz, and conveniently the suppressed carrier that is at 2441 MHz is in between channels.
These tones themselves can be modulated if necessary to spread the energy over the channel bandwidth, since single tone jammers can be easily defeated with tone excision. By spreading the bandwidth, it becomes a simple SNR competition and it can't be defeated given sufficient energy, as the duration of the signal occupation will be continuous in time.

If no sources were readily available, below shows a feasible implementation for comparison to alternative approaches. A custom 2 layer circuit board layout would be sufficient and is \$65 from Express PCB , (the sensitive RF areas should be minimized in pcb track lengths, and for this application there would be no issue implementing this design on such a pcb with attention to proper RF layout techniques in those areas)-  so a total cost of \$127 not including the low frequency noise source or DC supplies and DC control voltages. For single tones the low frequency noise source can be replaced with the same circuit as shown for the 51 MHz LO since the crystal oscillator (XO) can also be programmed for a 12 MHz output frequency. The noise source (modulated and filtered 12 MHz waveform) would be ideal for jamming purposes. The optional "carrier nulling" shown is a small DC offset (coupled in with a resistor) that can be adjusted to null any residual carrier feedthrough leakage at 2441 MHz if desired. Not shown in the spectrum above, but this would also result in suppressed carriers from the first upconverter stage in the output at 2413 MHz and 2467 MHz (also conveniently in between channels). The nulling control for this is in having a small DC offset control in the originating noise source. The nulling is feasible (and effective!) since the IF ports of both mixers have a response to DC.


Amplitude of first harmonic of a square wave, $2/\pi$ or $4/\pi$? I am getting really confused about the value of the first harmonic of a $50\%$ duty cycle $-1$ to $1$ square wave.
By doing the math I found $\frac{2}{\pi}$, in my lesson and Wikipedia it's $\frac{4}{\pi}$.
But now I checked with MATLAB and it's also $\frac{2}{\pi} \approx 0.6$
t = -100:0.1:99.9; % -> length = 2000
x = square(t, 50);
plot(abs(fft(x)/2000));

As you can see, on the spectrum the first harmonic is at $\frac{2}{\pi}$
And this is how I found $\frac{2}{\pi}$
\begin{align}
x(t) &= \text{rect}\left(\frac{t}{T/2}\right) \ast \delta_T(t) - \text{rect}\left(\frac{t}{T/2}\right) \ast \delta_T\left(t - \frac{T}{2}\right)\\
 &= \text{rect}\left(\frac{t}{T/2}\right) \ast \left( \delta_T(t) - \delta_T\left(t - \frac{T}{2}\right)\right)\\
\implies X(f) &= \frac{T}{2}\mathrm{sinc}\left(\pi f \frac{T}{2}\right)\cdot F \cdot \delta_F(f)\cdot\left(1 - \exp\left(-2\pi j f\frac{T}{2}\right)\right)\\
 &= \frac{1}{2}\mathrm{sinc}\left(\pi f \frac{T}{2}\right)\cdot\delta_F(f)\cdot\left(1 - \exp\left(-2\pi j f\frac{T}{2}\right)\right)\\
\implies X(F) &= X\left(\frac 1T\right) = \frac{2}{\pi}
\end{align}
Clarification and definitions I am using:
$$\text{FFT}\left\{\mathrm{rect}\left(\frac tT\right)\right\} = T\mathrm{sinc}\left(\pi f T\right)$$
$\delta_T$ is the dirac comb of period $T$.

It's just a question of the definition of the Fourier series coefficients. There is the complex Fourier series

$$x(t)=\sum_{n=-\infty}^{\infty}c_ne^{j2\pi nt/T}\tag{1}$$

where $c_n$ are the complex Fourier series coefficients, and $T$ is the period of the periodic function $x(t)$. For real-valued $x(t)$ you can also use the real-valued Fourier series

$$x(t)=a_0+\sum_{n=1}^{\infty}\left[a_n\cos(2\pi nt/T)+b_n\sin(2\pi nt/T)\right]\tag{2}$$

with

$$a_0=c_0,\quad a_n=2\text{Re}\{c_n\},\quad b_n=-2\text{Im}\{c_n\}\tag{3}$$

For a $T$-periodic square wave with duty cycle of $50\%$ and amplitude $1$ you get

$$c_n=\begin{cases}\frac{2}{jn\pi},\quad n\text{ odd}\\0,\quad n\text{ even}\end{cases}\tag{4}$$

From $(3)$ you get

$$a_n=0,\quad b_n=\begin{cases}\frac{4}{n\pi},\quad n\text{ odd}\\0,\quad n\text{ even}\end{cases}\tag{5}$$

So if you only use positive frequencies in your series respresentation, as in $(2)$, your get a factor of $4/\pi$ for the first harmonic, whereas if you use positive as well as negative frequencies, as in $(1)$, you get a factor of $2/j\pi$ for $n=1$ and for $n=-1$.

Your representation as a Fourier transform corresponds to $(1)$, where for each frequency component you get a Dirac impulse in the frequency domain. Since that Dirac comb has components at positive as well as at negative frequencies, the amplitudes equal $2/n\pi$:

$$X(f)=\sum_{n=-\infty}^{\infty}c_n\delta\left(f-\frac{n}{T}\right)\tag{6}$$


Denoising a signal using eigendecomposition I have a complex observable series $Y(t)$ which is the result of summing two complex r.v $X(t)$ (unobservable) and a $\epsilon(t)$ (unobservable).

$$Y(t)=X(t)+\epsilon(t)$$

Assume that $X$ and $\epsilon$ are independent, and that $\epsilon$ is serially uncorrelated and has constant mean and variance.

I construct a Hankel Matrix of the observed series and then I create a sort of variance-covariance matrix by multiplying the later by its conjugate transpose.

$$ R_Y=\mathcal{Y}\cdot \mathcal{Y^*}=\mathcal{X}\cdot \mathcal{X^*}+\mathcal{E}\cdot \mathcal{E^*}$$ (The last equality holds because of the independence of $X$ and $\epsilon$.)

After that, I make the eigendecompositon of that matrix and I obtain:

$$U\cdot (\Lambda_X +\sigma^2\cdot I)\cdot U^T$$  Where $U$ is the orthogonal matrix of eigenvectors, $\Lambda_X$ is a diagonal matrix containing the eigenvalues (in decreasing order) of $\mathcal{X}\cdot \mathcal{X^*}$ and $\sigma^2$ represents the "variance" of $\epsilon$ (under my assumptions $\mathcal{E}\cdot \mathcal{E^*}$ is a diagonal matrix with $\sigma^2$ in its diagonal.

Assume that the matrix $\mathcal{X}\cdot \mathcal{X^*}$ has not full rank, and because of that it has some null eigenvalues.

By the eigendecomposition of $R_Y$ we could select which eigenvalues correspond to the $\mathcal{X}\cdot \mathcal{X^*}$ matrix, (by considering the first $d$ eigenvalues, where $d$ represents the rank of the matrix $\chi \chi ^*$) and then we can estimate $\sigma^2$.

Knowing $\sigma^2$ we could estimate $R_X=\mathcal{X}\cdot \mathcal{X^*}$.

  Is there any way to factorize $R_X$ in order to obtain $\mathcal{X}$ and
  then by the  one-to-one relation between Hankel matrices and Series,
  obtain $X(t)$? if this is possible, under which additional conditions?

So I'll get this answer started by saying again that what you are describing reminds me a lot of Singular Spectrum Analysis.  The process is pretty much exactly how you describe: you have some observed signal, $Y[n]$, which is some sort of mixture, and we take $Y[n]$ and create a Hankel matrix from it using some desired frame length M (this will effectively determine your sub-space resolution).

Now that we have our Hankel matrix, we can either compute $R_Y$ directly and then do an eigendecomposition, or we could go with a Singular Value Decomposition of the trajectory matrix.  You'll find them to essentially be equivalent if we use the left singular vectors; either method is perfectly valid.

So once we get our matrix of eigenvalues, $U$, we need to project them using your typical Karhunen-Loeve method to get what we'll call principal components (this is sort of classic SSA notation, I dislike the naming but I digress):

$P = Y^H U$.

Now that we're here, we can use the singular values to determine amplitude of the signals themselves, and the principal components essentially determine the "parts" of the signal, i.e. $X[n]$ and $\epsilon[n]$.

Let's say we have N principal components, and we know that we only "need" the first one: from our principal components, we perform what is typically called Eigentriple grouping via the following multiplication:

$C = UP^H$

This new matrix, C, is what we call the reconstructive component(s) of $Y$, and is based on which principal components we selected.  We're not quite done yet, becuase $C$ is the same dimension as our Hankel matrix, $Y$, so do get back to our signal, we'll need an additional step.  To reconstruct the series itself, we'll perform what's called diagonal averaging (a bit of a misnomer I think).  The goal here is that we'll finally map back to a single signal, let's call it $Y_out[n]$.  It's in this step that we'll be able to utilize that one-to-one relationship of the Hankel to the series to extract out signals.  This step is a bit long winded, so I've omitted it for brevity, but you can consult Wikipedia or any of the many papers/texts on SSA for more information.

So now that we've briefly gone through it, you might find yourself asking "when will this method actually work".  SSA, like anything, is a tool and it has its limits.  From a signal processing standpoint, it's great at determining individual signals which are fairly narrowband (i.e. they don't cover the entire spectrum) and are "separable".  In the context of SSA, "separable" means that the signals are non-overlapping in frequency.  For a simple example, let's say I have two sinusoids, one at 50Hz, and another at 500Hz.  SSA will be able to easily separate the two signals, and we will be able to isolate the signals pretty easily; we would call these signals separable. However, let's say I change the second signal to be 51Hz instead of 500Hz.  In this case, SSA will struggle, as the signals are very, very close to one another; we would call these signals non-separable.

How do we determine which signals we care about? That's the difficulty with this method.  If you have some appreciable signal-to-noise ratio (SNR), you can perhaps use that to your advantage and identify which signals are of interest and which signals are just purely noise by inspecting the eigenspectrum; true signals as represented by the principal components will likely be associated with large eigenvalues, whereas noise will be relatively low power, and spread.

So, in short, if you have some time series observation, you can use SSA to attempt to discern individual signals/trends.  It's a really powerful tool, and as long as the signals are separable, you should have some reasonable success.  Hope that helps!


Why does my sinusoid look "AM" in shape? My code is :

Fs=200e6;
Ts=1/Fs;
NFFT=2^14;
Runtime=(NFFT-1)*Ts;
t=0:Ts:Runtime;
f_in=90*1e6;
y_in=sin(2*pi *f_in *t);
plot(t,y_in)
ylim([-1.5 1.5])

Then why does my plot look like amplitude modulated when you zoom into it?

This visual phenomenon appears because the maximum frequency is close to the Nyquist frequency, or half the sampling frequency. Sampling begins to approach the limit of $2$ samples per period, and thus the linear interpolation performed by Matlab becomes highly inaccurate. However, samples are correctly located, as you can see from the code where an higher sampling ('Oversampled') is superimposed:

Fs=200e6;
Ts=1/Fs;
NFFT=2^14;
Runtime=(NFFT-1)*Ts;
t=0:Ts:Runtime;
f_in=90*1e6;
y_in=sin(2*pi *f_in *t);

Fs2=20*Fs;
Ts2=1/Fs2;
NFFT=2^14;
t2=0:Ts2:Runtime;
f_in=90*1e6;
y_in2=sin(2*pi *f_in *t2);

clf;hold on;
plot(t,y_in,'x')
plot(t2,y_in2,'-')
ylim([-1.5 1.5])
xlim([5.2 5.4]*1e-6)
legend('Original','Oversampled')


Complex signal sine reconstruction Noob here, I read that any signal can be made by putting together sines and cosines, it always shows some kind of basic harmonic wave with constant amplitude such as square wave.

I understand that constant amplitude, constant frequency square and triangle waves can be made from sines, but what about square wave that sweeps its frequency or even suddenly jump to other frequency, what about square that has its amplitude changed?

What about noise, lets say I have 1 minute long 48KHz 16bit PCM signal that is just white noise, that can be reconstructed from bunch of sines too?

And what about transients, either the sudden square wave like ones or the gentle slow fade in type ones. Lets say I have signal that is silence and then square wave slowly and smoothly rises in amplitude, how can sines, which are constant in amplitude ever recreate it?

Basically, my point is these sines are constant in amplitude, frequency, and phase and they run the entire lenght of whatever signal we want to reconstruct, signals can have silence, transients, periods where amplitude and frequency is constant and periods where they change. I dont understand how can bunch of sines that run the entire lenght of complex signal ever reconstruct it. How can 1000 sine waves sum up into perfect silence and then suddenly sum into noisy square wave sweep?

A sum of sines/cosines (with frequencies that are integer multiples of a fundamental frequency) will always result in a periodic function (that's what Fourier series are about). But you can of course approximate any function with finite support by choosing the length of the function's support as your period (if you don't mind that the approximation will be a periodic continuation of the given function).

Note that the approximation by a Fourier series is a least squares approximation, not a point-wise approximation. So if the desired function has discontinuities there will always remain some error of fixed magnitude, no matter how many sinusoids you use for the approximation. This oscillatory behavior near discontinuities is called Gibbs phenomenon.


Preparing audio data for FFT I would like to experiment with some input from a microphone and am receiving a 2 channel, 512 samples buffer in real time

I know the signal could be passed through a low pass filter, and then windowed before the FFT.

The sample rate is 44100.00, what low pass filter is needed for this, does the analog to digital conversion not run it through a low pass filter before converting it?

How exactly does a window apply to real time sample input of 512 samples? Is a window not the same as a filter?

I have also read that the signal can be run through the FFT twice, why exactly is this done and is it necessary for this data?

I would like to slowly build this up to do spectral analysis etc. Is this done directly on the FFT bin data?

I know the signal could be passed through a low pass filter, and then windowed before the FFT.

Yes, it could be, and it's usually advisable to use a window, but these are not just things you do blindly because you heard about them. Low passing the data will alter it in ways you might not want to. It might be useful, or it might not -- depends on your application. It will reduce the amount of high frequencies represented in your FFT results. Is that what you want? Maybe, maybe not. Windowing is useful because it reduces (but does not eliminate) frequency-domain distortion caused by separating your data into chunks. Usually you want that, but there are trade-offs (my explanation is an over-simplification).

  The sample rate is 44100.00, what low pass filter is needed for this, does the analog to digital conversion not run it through a low pass filter before converting it?

The input must be filtered before the digital conversion (and usually this is done for you). There is no way to filter it after to achieve the same effect, so even if it's not taken care of for you, you might as well forget about it. Chances are, your DA Converter does a halfway decent job of this, though.

  How exactly does a window apply to real time sample input of 512 samples? Is a window not the same as a filter?

A window and a filter are different, as explained above.

  I have also read that the signal can be run through the FFT twice, why exactly is this done and is it necessary for this data?

You haven't really told us anything about what you want to do with this data, but no, it's highly unusual to do two FFTs.

  I would like to slowly build this up to do spectral analysis etc. Is this done directly on the FFT bin data?

More or less. I suggest you read up on filtering, windowing, fourier analysis on Wikipedia. This blog post might also help.


How to filter pulse like jumps in a signal? I have a noisy signal with jumps to high value and back to normal for limited time bins (40 for example) which I want to filter to have a smooth signal with no jumps. Here is the graph, the blue one.

I tried moving average filter (n=20), the green one. Also tried window based (kaiser) LPF with cut-off freq=0.01 (I choose it manually) but still the filtered signal follows the jumps.

Also I tried to flatten the signal by using the average of the past n values whenever there is a jump. But how could I find the start of a jump. I need a threshold. Also there is difficulty with determining the end of a jump. I need another threshold.

Is there any way to filter the signal in a way that it doesn't follow
the jumps?
Are those jumps and falls caused by low frequencies, since LPF don't eliminate them? So I should use high pass filter? I am confused.

EDIT: Another graph that shows what I want to eliminate (the gray regions).

I think that you will not be able to get acceptable results with linear filtering because the impulse noise covers the whole frequency range. You may want to look into nonlinear filtering techniques. As a first (and simple) approach try median filtering: http://en.wikipedia.org/wiki/Median_filter

I think if you optimize the window size of the median filter you should get better results than with purely linear filtering.


How are two signals multiplied? And how is it different from convolving two signals? Suppose that I have two signals $x[n] = \left\{2,4,1\right\}$ and $p[n] = \left\{5,1,8\right\}$ and I want to multiply them.

How do you do that?
How different is it from convolving two signals?

I understand that multiplication in one domain is equal to convolution in other domain. How do you choose as to what to use : multiplication or convolution?

Yes, you are correct. Multiplication in time domain means convolution in frequency domain and vice versa.  Multiplying your signals $x[n]$ and $y[n]$ will give an output:

\begin{align}
z[n]&=\{2\cdot 5, 4\cdot 1, 1\cdot 8\}\\ &= \{10, 4, 8\}\end{align}

Remember that this output is in time domain. When you convolve $x[n]$ and $y[n]$, you will get $z[n]$ in time domain as:

\begin{align}
z[n]&=\{5\cdot 2, 5\cdot 4+1\cdot 2, 5\cdot 1+1\cdot 4+8\cdot 2, 1\cdot 1+8\cdot 4, 8\cdot 1\}\\ &= \{10,22,25,33,8\}
\end{align}

Just flip one of the signals around zero and start moving right one place at a time. Multiply the corresponding points as you go along. The output has a larger sequence because convolution output has ${\rm length}(x)+{\rm length}(y)-1$ points.

To answer your query about where to use multiplication & convolution, assume you want to pass signal $x(n)$ through filter $y(n)$. The output of the filter $z(n)$ will be convolution of $x(n)$ and $y(n)$.
Now assume that you first converted from time to frequency domain, i.e. $X(e^{i\omega})$ and $Y(e^{i\omega})$ are frequency domain representation of $x(n)$ and $y(n)$. Now, to find $Z(e^{i\omega})$(output in frequency domain), you have to multiply $X(e^{i\omega})$ and $Y(e^{i\omega})$. To get the output in time domain i.e. $z(n)$ you have to apply inverse transform.

So that's how you use convolution and multiplication.


Is the following system stable and causal? Suppose the following function describes the unit step response of a system, where $u[n]$ is the unit step function.
$$
y[n]=\left(\frac{1}{2}\right)^{n-1}u[n+1]
$$
I want to find out the system function and if the system is stable and causal. My first step is to calculate $Y(z)$.
$$
\begin{align}
Y(z)&=\sum_{n=-\infty}^{\infty}y[n]z^{-n}\\
&=\sum_{n=-\infty}^{\infty} \left(\frac{1}{2}\right)^{n-1}u[n+1] z^{-n}\\
&=2\sum_{n=-1}^{\infty}\left(\frac{0.5}{z}\right)^n\\
&=2\left( \frac{1}{1-0.5z^{-1}}+\frac{z}{0.5} \right)\\
&=\frac{4z^2}{z-0.5}
\end{align}
$$

Above I exploited the fact that I can form a geometric series which converges for $|z|> \frac{1}{2}$. Now calculating $H(z)$ is straight forward, since we know the input is $u[n]$.
$$
H(z)=\frac{Y(z)}{X(z)}=\frac{\displaystyle \frac{4z^2}{z-0.5}}{\displaystyle \frac{z}{z-1}}=\frac{4z(z-1)}{z-0.5}
$$

We see that the system has two zeroes,  $z_{0,1}=0$ and $z_{0,2}=1$, as well as one pole at $z_{\infty,1}=0.5$.

My last step was calculating the impulse response for this system.
$$
\begin{align}
H(z)&=\frac{4z^2}{z-0.5}-4\frac{z}{z-0.5} \\
h[n]&=\left(\frac{1}{2}\right)^{n-1}u[n+1]-4\left(\frac{1}{2}\right)^nu[n]
\end{align}
$$

From the above calculations I want to find out if the system is stable and causal. To see if it is causal I need to check if the impulse response disappears for negative time indices. Well, that is not the case, since $h[-1]=0.5^{-2}=4$. That means we have an acausal system, where the region of convergence (ROC) lies in between $|z|=0$ and $|z|=\frac{1}{2}$. Therefore, the system is also not stable, since the unit circle is not in the ROC.

Here I am getting confused: The system converges for $|z|<\frac{1}{2}$, but for calculating $Y(z)$ $|z|$ must be greater than $\frac{1}{2}$. How can that be? Are my calculations and thought process correct?

You are on the right track, but you are confused.
First, regarding the calculations you could simply do
$$y[n]=2\left(\frac{1}{2}\right)^n(u[n]+\delta[n+1])=2\left(\frac{1}{2}\right)^nu[n]+4\delta[n+1]$$
which gives you

$$Y(z)=\frac{2z}{z-\frac{1}{2}}+4z,\,|z|>\frac{1}{2}$$
The impuse resposne is the first difference of the step response:

$$\begin{align}
h[n]&=y[n]-y[n-1]\\
&=2\left(\frac{1}{2}\right)^nu[n]+4\delta[n+1]-4\left(\frac{1}{2}\right)^nu[n]\\
&=4\delta[n+1]-2\left(\frac{1}{2}\right)^nu[n]
\end{align}$$
Obviously, $h[-1]=4$ and the system is not causal.
But here you are mixing up finding the inverse $z$-transform of anti-casual signals with system being non-causal. ROC is $|z|>\frac{1}{2}$ and the system is non-causal full stop.


Understanding overlapping in STFT I understand the concept of the STFT. In order to avoid spectral leakage, you use a Hann window that overlaps by 50%. However, the problem that I have, is that the example code online does not account for this.

What I mean is, if the signal is overlapping by 50%, why doesn't anyone add the frequency content of the overlapping sections?

https://kevinsprojects.wordpress.com/2014/12/13/short-time-fourier-transform-using-python-and-numpy/

This is an example of what I mean. It looks like he is grabbing the overlap, but he is not adding the overlap pieces together. I don't see how the frequency content is being saved. He windows a function, does the FFT, then moves over and does the FFT again with the overlapping piece. I feel like he should be adding the overlapping pieces together. Is there something I am missing?

Thank you: John

Edit: I see the first answer, but it is still not hitting home. Hanning windows make the left and right information approach zero. So when it is overlapped, this removes this problem. However, when the FFT is taken, the information of the transform is put side-by-side. Wouldn't this mean, First Window=> upper time series information lower amplitude, Second window =>lower time series information lower amplitude.... This does not seem to solve the issue. You are still losing information with each window... is this not the case?

I understand the concept of the STFT. In order to avoid spectral leakage, you use a hann window that overlaps by 50%.

I'm sorry but you have a misunderstanding of spectral leakage in addition to how a spectogram should be computed. To be exact you cannot avoid spectral leakage completely; all you can do is to make a compromise between the spectral resolution and the spectral leakage: A window that yields less spectral leakage should have smaller side lobes (in its window Fourier transform), but that window will unavoidably smear the input signal spectrum because of the increased main lobe width (of the window Fourier transform) and therefore limit the spectral resolution per given window length.

So you should optimize between spectral leakage and spectral resolution. In order to increase spectral resolution for a given accepted leakage level you can increase the window length. But then this will decrease the time localization of the STFT. And in order to accomadate for this latter side effect, you can use overlapping windows so that the next window does not start from the end of the current window, but begins somewhere in the middle; an accepted choice is 50% overlapping. But you can also use other percentages according to your application as @StanleyPawlukiewicz states in his answer.

Note that 50% overlapp is also an optimal choice to reduce the variance of the estimate of the power density spectrum computed by the periodogram averaging applied in overlapping windowed analysis segments (the Welch's method)


Relation between samplingrate and frequency I am working on Fourier Transformation, and applying this for recognizing an audioclip.

I have a 9 second long audio clip of a guitar strumming an A-Minor.
The audioclip has a sampling rate of 44100 Amplitudes per second, and a bitdepth of 16 bits.

I do a fourier Transform of the audiofile(green graph), and get a spectrum where the samples are the x-axis. I want to have frequencies along the x-axis, for being able to recognize the pitches, so I again can recognize the chord that is being played.

Instead I get samples. How can I convert this to frequency?

Here is a photo of my spectrums, one of them zoomed in on the samples from 0 to 10 000.

you have 9 $\times$ 44100 = 396900 samples. I assume you did a single FFT and you should take the first , 198450 complex points for your plot after taking the magnitude square.  The last 198450 points out of the FFT are a mirror image of the first 198450 points so you can neglect them. The first kept  point corresponds to $f=0$ and the last kept point to $f=22050$ Hz.  which makes each bin  .1111 wide


Confusion in Normalised Bandwidth of a signal I know normalised bandwidth can be found from FFT plot of the signal. But I'm confused that at what stage I have to take the FFT.

I mean first binary data is input to a convolutional encoder, then to NRZ encoder, Then GMSK modulation (sampling frequency till now is 100KHz), then upsampling to 30MHz sampling frequency and then the quadrature modulation.

At what stage I'm supposed to take the FFT to find normalised bandwidth?

You can take the FFT on any signal after modulation, but would be simplest (fewer samples required) to do it immediately after modulation. With proper upsampling the FFT of your signal prior to upsampling would be identical as long as you used proportionally more points to keep the time length identical. In your example of 100 KHz and 30 MHz, you could take a N point FFT at 100 KHz or a $300 N$ point FFT at 30 MHz and the result should look identical (with just more points interpolated in the FFT due to the upsampling). Given the 100 KHz sampling rate and room for your interpolation filter, I would assume that your GMSK modulated signal BW you are trying to view is less than 70 KHz (as a complex signal)?


Notation of an LTI system consisting of LTI filters I would like to find a reference for two notations of an LTI system consisting of LTI filters. In z-domain, the LTI system is given by

$$
 \mathbf{y}(z) = \mathbf{C}(z) \mathbf{s}(z) + \mathbf{D}(z) \mathbf{x}(z) \\
 \mathbf{s}(z) = \mathbf{A}(z) \mathbf{s}(z) + \mathbf{B}(z) \mathbf{x}(z) ,
$$
where $\mathbf{x}$ and $\mathbf{y}$ are the input and output of the system, respectively and $\mathbf{s}$ is the system state. The matrices $\mathbf{A}$, $\mathbf{B}$, $\mathbf{C}$ and $\mathbf{D}$ constist of LTI filters, e.g.,
$$
A_{ij}(z) = \frac{b_0 + b_1 z^{-1} + \dots + b_k z^{-k}}{a_0 + a_1 z^{-1} + \dots + a_k z^{-k}}.
$$
Now, each LTI filter can be represented by its impulse response in time domain. Let the matrices $\mathbf{A'}(n)$, $\mathbf{B'}(n)$, $\mathbf{C'}(n)$ and $\mathbf{D'}(n)$ with $n$ indicating time such that e.g.
$$
A'_{ij} = [a'_0, a'_1, a'_2, \dots ]
$$
is the impulse response of the LTI filter $A_{ij}(z)$.

Is it then possible to rewrite the whole LTI system in time-domain as
$$
 \mathbf{y}(n) = (\mathbf{C} \ast \mathbf{s})(n) + (\mathbf{D} \ast \mathbf{x})(n) \\
 \mathbf{s}(n) = (\mathbf{A} \ast \mathbf{s})(n) + (\mathbf{B} \ast \mathbf{x})(n) ,
$$
where $\ast$ denotes the convolution operation. I would especially appreciate any reference to literature / text books where this representation is discussed.

An one-dimensional example:
The state-transition matrix may be a (time-invariant) one-pole filter, e.g., $\mathbf{A}(z) = \frac{b_0}{1 - a_1 z^{-1}}$. The corresponding time-domain $\mathbf{A}'(n)$ is the impulse response of this one-pole filter.

As I've said in chat, I do not believe the system
$$
 \mathbf{y}(z) = \mathbf{C}(z) \mathbf{s}(z) + \mathbf{D}(z) \mathbf{x}(z) \\
 \mathbf{s}(z) = \mathbf{A}(z) \mathbf{s}(z) + \mathbf{B}(z) \mathbf{x}(z) ,
$$
with the matrices time-varying is an LTI system.

That means, I don't think we can answer your question.

If the matrices are not functions of $z$ then you can write the transfer function as:
$$
\mathbf{y}(z)/\mathbf{x}(z) = \mathbf{C}(zI - \mathbf{A})^{-1} \mathbf{B} + \mathbf{D}
$$
but it's not clear to me if this equation is valid if the matrices are time-varying / $z$-transforms.


difference in spectrum First I defined 5 signals with same length (40 samples) with different frequency and sampling frequencies, then I've put them together like [X1 X2...] and took it's fourier transform and the result is in the corresponding figure, the I took those 5 signals with same length (200 samples) and sum the up and took it's fourier transform, the result is in the corresponding figure, what makes these two signals fourier transforms different? why is it?

If I understand correctly, the OP is comparing cascading 5 time domain signals one after the other as in $[f_1, f_2, f_3, f_4, f_5]$ and taking the DFT to summing longer samples of the signals and then taking the DFT. The following is what we would expect to see in this case.

The DFT of the sum of functions is equal to the sum of their DFT's:

$$\sum_{n=0}^{N-1}f_1(n)W_N^{nk}+\sum_{n=0}^{N-1}f_2(n)W_N^{nk} = \sum_{n=0}^{N-1}(f_1(n)+f_2(n))W_N^{nk}$$

Consider first all sequences of the same length and cascading versus summing the sequences:

Cascading $R$ sequences of equal length $M$ is identical to zero-padding each of the sequences out to $N = RM$, each with  starting offset of $rM$ where $r$ is from $0$ to $R$.

Example with two sequences with $M=4$:

$f_1 = [1, 2, 3, 4]$

$f_2 = [6, 8, 9, 1]$

$\text{DFT}\{[1, 2, 3, 4, 6, 8, 9, 1]\}$

$ =  \text{DFT}\{[1, 2, 3, 4, 0, 0, 0, 0]\} + \text{DFT}\{0,0,0,0, 6, 8, 9, 1]\}$

Zero-padding a sequence prior to taking the DFT does not change the original DFT values, but inserts new interpolated values in between the samples in frequency. Such that

$$F_p(nk) = F(k) \tag{1}\label{1}$$

Where $F$ is the DFT of a sequence $f$ and $F_p$ is the DFT of $f$ with zero-padding. This is proven below along with the fact that equation $\ref{1}$ holds true for all offsets $rM$.

So in the example with $f_1$ aand $f_2$ above, $$F[k] = \text{DFT}\{[1, 2, 3, 4, 6, 8, 9, 1]\} $$

$$F[nk] = \text{DFT}\{[1, 2, 3, 4]\} + \text{DFT}\{[6, 8, 9, 1]\}$$

Similarly, the DFT of the cascade of the five sequences would have the same values as the DFT of the sum of the five sequences with four interpolated values in between each of the DFT samples!

Since in this case the OP didn't simply sum the sequences, but rather summed longer sequences out to the same length as the cascaded sequence, then here the DFT of the longer sequences would have the higher frequency resolution and scaled 5 times larger given by the longer data set, so the values in between that were previously interpolated would now be dictated by the longer time samples as new DFT values. However every 5th sample in both sets should match exactly other than the scale of 5 for the same reason outlined above, with more samples of actual frequency data in between!

Proof of relationship used above in equation $\ref{1}$

Given sequence $f(n)$ where $f(n) = 0 \text{ for } n \notin [0, M-1]$

The DFT of $f(n)$ is:
$$F(k)=\sum_{m=0}^{M-1} f(m)W_M^{mk}, \space\space\space\space k = 0,...M-1\tag{2} \label{2}$$

Where $W_M^{n}$ are the $n$ roots of unity given by $e^{-j2\pi n/N}$

The zero padded DFT of $f(n)$  out to length $N=RM$ with $R\ge 0$ with any offset $rM$ with $r \in [0, R]$ is:

$$F_P(k) = \sum_{n=0}^{N-1} f(n-rM)W_{N}^{nk}, \space\space\space\space k = 0,...N-1\tag{3} \label{3}$$

Using the DFT shift theorem, equation $\ref{3}$ becomes:

$$F_P(k) = \sum_{n=0}^{RM-1} f(n-rM)W_{RM}^{nk}=W_{RM}^{rMk}\sum_{n=0}^{RM-1} f(n)W_{RM}^{nk} $$
$$= W_{R}^{rk}\sum_{n=0}^{RM-1} f(n)W_{RM}^{nk}$$

Since $f(n) = 0$ for $n>M-1$

$$F_P(k) = W_{R}^{rk}\sum_{n=0}^{M-1} f(n)W_{RM}^{nk}\tag{4} \label{4}$$

$$F_P(Rk) = W_{R}^{rRk}\sum_{n=0}^{M-1} f(n)W_{RM}^{nRk} = \sum_{n=0}^{M-1} f(n)W_{M}^{nk} = F(k)$$

Thus we see how given any sequence $f(n)$ of length $M$ with DFT $F(k)$, zero padded out to length $N=RM$ for any positive integer $R$ with DFT $F_P(k)$, the following equality will hold for any rotational offset (starting index) of $rM$ for any integer $r$ from $0$ to $R-1$:

$$F(k) = F_p(Rk) \tag{5} \label {5}$$


Gesture recognition with Hidden Markov Models I am currently working on a gesture recognition system (for an Android Application). I think that I have completed the Image processing stage, where I am able to extract the contour of the hand (I am wearing a glove to avoid background subtraction for now).

I am also getting the bounding ellipse/rectangle, the centroid as 'important' shape features of the hand.

My problem is that not a lot of literature indicate what the next stage is before the classification of the image through the hidden markov model. I am finding this 'Feature Extraction' stage very ambiguous.

Currently I am getting a list of angles of the contour (which is approximated so as to obtain a limited number of angles)

The problem is that I am clueless as what to do next. When I search for literature as 'Shape classification using HMMs' I still find it hard to what should I do next.

As a tool for HMM I plan to use JaHMM. But I don't know how to experiment with the tool at this stage since I do not know what will be the input to this classification stage!

I have encountered a list of functions I found in some literature, such as Fast Fourier Descriptors, Curvature Descriptors, B Spline; yet I am clueless as to how to apply these functions to my current data (i.e. list of angles, e.g.: -63, 154, 3, 23, 54, ....)

UPDATE 1:

Thank you for your information. @Peter K.

Regarding the poses: I was going to go for a set of words in American Sign Language which are fingerspelled, e.g 'dog' -> 'd' 'o' 'g' (3 states); but the moment I haven't decided what to go for.

I am currently reading some more papers to see what type of information I should extract, such as :

centroid of hand
angle of motion
distance from a particular point to different points of the hand contour (etc..)

Now I have encountered a paper which seems to show what I want to do, I'm not sure:

http://espace.library.uq.edu.au/eserv.php?pid=UQ:10700&dsID=n0273.pdf

I am currently reading section 5 - Vector Quantization (I heard of this term before but do not actually know what it involves, and in figure 5 there seems to be a complex algorithm which, if I understood correcty, converts the set of values I achieve from the hand (just mentioned) into 1 digit which I can use as the Observation sequence to train an HMM for 1 particular sign. Do you think I'm moving on the right track? (I'm working on Android, (NDK), I found JaHMM as an HMM tool, and using OpenCV for image processing.

@Peter K. Thank you for your answer. Regarding the data generation I was planning to follow the steps of this paper, after I produce my personal dataset:
(section 4 and 5)
http://www.i.ci.ritsumei.ac.jp/~shimada/papers/vi02-tanibata.pdf

UPDATE 2: keeping in mind that a gesture is formed of {posture a, posture b, posture c}

I am now thinking that I must make use of some form of classification algorithm. That is, I currently have a set of feature vectors:

Posture A:
[angle of ellipse surrounding it, height:width ratio]

0.802985 33.909615
0.722824    31.209663
0.734535    30.206722
0.68397 31.838253
0.713706    34.29641
0.688798    30.603661
0.721395    34.880161

Posture B:
[structured the same as posture A]

0.474164 16.077467

0.483104 14.526289
0.478904 14.800572
0.483134 14.523611
0.480608 14.41159
0.481552  15.563665
0.497951  15.563585

etc..

and I would like that when I feed a feature vector I obtain a simple symbol, e.g. 'A', 'B', etc.

Is this possible? I also migrated the question here: https://stackoverflow.com/questions/15602963/vector-quantization-algorithms-used-to-provide-observation-sequences-for-hidden

I have used HMM for gesture recognition (not pose recognition). What I did was: tracking the hand and recognize the gesture the hand was drawing in the air, you can image it as a trail.

You can use HMM as sequence recognizer, so first of all you need to transform your image into a discrete number sequence.

For each gesture you want to recognize, you need to train an HMM for that.

So you have a dictionary with some known word. Each one is a trained HMM. If you have a new word (unknown obdervation) you can compute the probability for each word of the dictionary to be likely the unknown one.

Some pseudo-code:

##################### training phase

examples = [112233, 11233, 123, 1122223333]
word1 = train(example)

examples = [222333, 22222223333, 2222333333]
word2 = train(example)

examples = [124555, 1122445, 1111111222224444555]
word3 = train(example)

dictionary = [word1, word2, word3]

##################### recognition phase

#let's say you have a new unkown word: 12245, you want to know what word of the dictionary it is more likly

unkown = 12245
probabilities = []
for w in word:
    probabilities.append( calculate_prob(unkown) )

Now, see what the is the maximum value in probabilities and you get the most likely word of the dictionary!

See here:

http://www.ece.ucsb.edu/Faculty/Rabiner/ece259/Reprints/tutorial%20on%20hmm%20and%20applications.pdf

http://www.codeproject.com/Articles/69647/Hidden-Markov-Models-in-C

http://www.creativedistraction.com/demos/gesture-recognition-kinect-with-hidden-markov-models-hmms/


Measuring the average frequency of signals I have many signals/audio data of a hammer banging on a nail. As the nail gets deeper it seems obvious to me (just by hearing) that the pitch is higher at the end than at the beginning. I assume that the pitch gets higher with every bang. I want to be able to measure this.

So my question is how to I compute the frequency or average frequency of each bang, to be able to see this increase in pitch?

What I thought about doing was take the fourier transform of each bang and then take a weighted average. Like: (sum over (fourier coefficient index * amplitude of that coefficient)) / number of fourier coefficients.

But I am not seeing a clear rise in this value. Perhaps what I should care about just the fourier coefficeint with the largest amplitude?

What methods (if mine is not correct) to get the average frequency?

The weighted average "(sum over (fourier coefficient index * amplitude of that coefficient)) / number of fourier coefficients" does not seem very good as it will always give more weightage to the higher frequencies( because fourier coefficient index for higher frequencies will be higher) and the frequencies which you are looking for might be having low index.

Better method could be to see which fourier coefficients have maximum values for each bang of hammer and corresponding to that you can find the major frequencies(band of frequency) for each bang.

hope this helps.


What Are Different Approaches to Realize a Gaussian Blur (Smoothing) Step on an Image? Could some review some methods to apply Gaussian Filter (Blur) on an image besides the direct one using Truncated FIR (classic convolution with Truncated Kernel) approximation?

You can apply Gaussian Blurring on an image in many ways:

Using FIR Approximation by Convolution.
Using Approximation by Box Blur.
In the Fourier Domain by Multiplication by a Fourier Kernel.
Using IIR Approximation.

You may have a look on my project - Fast Gaussian Blur.

Update

Meaning of FIR

In this context FIR is the coefficients which are utilized to create the convolution kernel.
Indeed in the case of Gaussian Blur the 2D Kernel can be represented by outer product of 2 1D kernels hence it can be done in a separable fashion - Rows / Columns and then Columns / Rows.

Multiplication in The Fourier Transform

By utilizing the Convolution Theorem which states there is equivalence between Convolution in Spatial Domain to Multiplication in Fourier Domain one could use it to apply the convolution.

Run Time

This really depends on 3 things:

Size of the image.
Size of the kernel (The STD of the Gaussian Kernel).
The Quality required.

For instance for cases of large image and large kernel (Size which is approximately the size of the image) the Fourier Transform is really good especially in the case we have few images to be smoothed with the same kernel (For instance 3 RGB Channels).

Usually, in practice, in most times, the methods which their complexity is independent the Kernel STD (Box Blur Approximation / IIR Approximation) are the fastest methods with acceptable quality.
You may try them in the project I linked to above.


Fourier transform of dirac comb with function: The scaling factor Multiplication in the time domain corresponds to convolution in the frequency domain:

$$
f(t) \cdot x(t) \iff F(j \omega) * X( j \omega)  \tag*{No scaling factor}
$$
I know the fourier transform of the dirac comb is:
$$
\mathcal{F} \big \{ \text{III}_{T_{s}} (t) \big \} =  \omega_{s} \cdot \text{III}_{T_{s} } (j \omega)
$$
But according to oppenheim and others:

$$
\mathcal{F} \big \{ f(t) \cdot \text{III}_{T_{s}} (t) \} =  \underbrace{\dfrac{1}{T_{s}}}_{\text{scaling factor}}  \cdot \displaystyle \sum_{k \to - \infty}^{ k \to \infty} X(j( \omega - k \omega_{s} ))  \text{ (Scaling factor)}
$$
So my question is, for the above, why does the scaling factor not remain $$ \omega_{s} $$ why does it become $$ \dfrac{1}{T_{s}} $$
Which means that
$$
x(t) \cdot y(t) \iff \dfrac{1}{2 \pi}  X(j \omega) * Y(j \omega) \tag{ ??? }
$$
Why does the scaling factor not remain the sampling angular velocity?

Your first formula is wrong. If you use angular frequency ($\omega$), then multiplication in the time domain corresponds to $1/(2\pi)$ times convolution in the frequency domain, just as in the last formula in your question:

$$\mathcal{F}\{x(t)\cdot y(t)\}=\frac{1}{2\pi}X(j\omega)\star Y(j\omega)\tag{1}$$

If you use frequency $f=\omega/(2\pi)$, then you get

$$\mathcal{F}\{u(t)\cdot v(t)\}=U(f)\star V(f)\tag{2}$$

where I've used the definitions

$$X(j\omega)=\int_{-\infty}^{\infty}x(t)e^{-j\omega t}dt\tag{3}$$

$$U(f)=\int_{-\infty}^{\infty}u(t)e^{-j2\pi f t}dt\tag{4}$$


What is the point spread function and optical transfer function and what uses are they in image processing So I seem to really be struggling with the concepts of the point spread function (PSF) and optical transfer function (OTF) and what they are used for in image processing.

Everything I seem to google on the concepts seems to come back with telescopes and microscopes or I seem to get information about Delta functions and impulse response of a system and how this relates to convolution, but nothing seems to tell me about what PSF or OTF actually do or measure.

From what I have gathered so far is the following.

You start with a linear time invariant (LTI) system you then apply a signal which can be broken down into Diracdelta functions ($\delta(x)$), aka the impulse. The (LTI) then outputs a response aka the impulse response ($h(x)$) and according to some web sites this is the PSF. Now I can follow the maths to a degree and see how the convolution between the impulse and the impulse response, but what I can wrap my head around is what the PSF is actually useful for.
Why is it useful for imaging systems? From a video that I glimpsed at on you tube the PSF spreads the Dirac delta function so in terms of image processing is this effectively telling us how much blur will be in the image itself? Is the the only thing the PSF tells us?
Moving on to the OTF, so this is the Fourier transform of the PSF and that as far as I have gotten, as I don't fully understand what the PSF dose for an imaging system. I really cant figure what the OTF does for an imaging system.
What's more confusing is now, I keep thinking of the OTF as a kinda of filter function when I see the convolution equation in the frequency domain, dose the OTF itself play a role in frequency filtering if so how?

Fundamentally could someone explain the concepts of PSF and OTF and what they actually tell me about the system in terms of image processing.

PSF (point spread function) is the spatial domain impulse response of the optical imaging system assuming it's modeled an LSI (linear shift-space invariant) system according to its light intensity processing.

In principle, PSF can obtained by inputting an impulse (varry small bright dot - a point) to the system, and measuring the output response where the dot is typically spread into a blurry circle; hence the name point-spread function.

Also for CRT display systems, PSF can be obtained by sending a single pixel dot signal to the display, and observing its resulting brightness spread on the screen.

Once you obtain the impulse response (or PSF) of an LSI system, then you can mathematically compute the output response of the system to any input by convolving the PSF with the input signal. This concept lies at the heart of the signal processing, and that's why PSF is a very important characterisation of optical systems.

OTF is the Fourier transform of the PSF, and therefore carries the same information, but, stated in frequency-domain instead.


Sub nyquist sampling, required number of samples for time sparse grouped signals Question: Does it make sense to perform compressed sampling if the non zero samples are grouped in time? If so, what is the minimal length of the vector x that should be acquired to allow full signal reconstruction (of at least the non zero pulse segments)?

$y=\phi x$, where y is the compressed vector and x is the signal (sampled at nyquist)

Situation: time sparse signal (+- 97 % of the samples at noise level), but the 3% non zero samples are grouped (periodic). I'm only interested in the pulse samples (so not in the noise floor samples).

Thanks a lot!

Compressed sensing does not assume any distribution of non-zero elements in the input vector (signal), so it makes no difference if your non-zero elements are near each other or located uniformly on the sampling grid. The number of measurement, as theory states for all sparse signals, depends on how many non-zero elements you have (sparsity degree), what is your noise in input and what recovery quality your need (exact or approximation) and etc.
therefore,

Does it make sense to perform compressed sampling if the non zero samples are grouped in time?

Yes!
I suggest the following video lecture of Baranuik:
Video Lecture On CS


Smoothed Square Wave new here and probably pretty unexperienced compared to the rest of you. This should be simple enough but just wanted some clarification. I'm trying to model an analog square wave oscillator with C. With the oscillator I'm trying to mimic, the square wave isn't totally square and has curves with the starting edge of each pulse (on/off). It looks like in the analog device the capacitor is smoothing the square wave, so I believe I need to implement a low pass filter in my code to smooth out the square wave.

tldr - Just wanted to confirm that I should implement a low pass filter on the raw square wave to make it sound similar to the analog oscillator of the picture below. If there's another way I should be approaching this I'm open to listen to suggestions as I might be approaching this wrong.

Top is my square wave
Bottom is device's square wave

Yes. This looks like it's a typical first order low pass. The time constant can be determined by looking at the time it takes for the falling edge to drop to 37% of the max amplitude ($e^{-1}$). The cutoff frequency of the low pass filter is $1/(2 \cdot \pi \cdot t_c)$ where $t_c$ is the time constant


Is sound the best way to transmit a signal under water? Could sound be used in transmitting sound waves under water digitally? Could the sound then be converted by a repeater to transit by radio waves? What other ways can a signal be transmitted other than sound over long distances under water?

https://space.stackexchange.com/questions/32421/how-would-a-transmission-be-relayed-through-bottom-of-europa-ocean-to-earth

Yes, it is called acoustic communications. Iterative Carrier Frequency Offset and Channel Estimation for Underwater Acoustic OFDM Systems is an example of a paper that uses orthogonal frequency division multiplexing (OFDM) in an underwater acoustic channel.
EDIT:
Note that you wouldn't call it a SONAR any more because SONAR stands for SOund Navigation And Ranging, whereas this is a communication system, just like you wouldn't call your cell phone a radar.
EDIT 2:
Yes, the signal could be converted to a radio signal at a repeater. The repeater would ideally be floating in the surface and have a transducer (which is basically a microphone or array of microphones) on the underwater side to receive the sound waves. The transducer converts the sound waves into an electrical signal, which can be filtered and/or amplified, and then sent to an antenna for radio transmission. The antenna would likely need to be above water because radio waves don't propagate very far through water.


Audio transfer function modelling with pink noise + signal I'm trying to model a transfer function for a noisy audio system, specifically to measure delayed system response. Before I can confidently apply control, I need to verify that I can exert control over the system by broadcasting a signal and listening for the delayed (and potentially noisy) feedback.

I suspect I can add or subtract a unique signal to pink noise to simultaneously calibrate my transfer function as well as verify that the system is responding to my input.

Is there any prior art on this? Is pink noise the right approach here?

Current Plan:

Generate a second of pink noise and calculate its FFT.
Generate a second of unique signal and calculate its FFT. Maybe this can be a simple square wave at 5kHz.
Sum the noise and carrier signal together as test signal. Split the signal into five 200ms chunks, calculate the FFT of the noise
Record a second of (noisy) silence.
Start broadcasting the combined test signal, while still recording.
Take the FFTs of the "empty" and "full" signals. The difference should be the ambient noise.
Look for the source spectrum distribution.
Given that we found the distribution we were looking for, cross-correlate the recording with the source signal to find the delay at k samples.

Does this sound about right? When I'm done, I should know the average ambient noise, the system delay and the harmonic distortion (maybe I need to sweep frequency for this?).

For calibration, performance doesn't have to be real-time. Is there a more complicated sweep function I can use to find a more accurate delay and response?

your generated pink noise is just a signal.  it's not noise because you know what it is.

you can toss whatever broadbanded signal you want into your linear, time-invariant system.  call that "$x[n]$".
you can measure the output of that LTI system, call that "$y[n]$", given the known input.
you can FFT both the input and the output.  call those FFT results "$X[k]$" and "$Y[k]$", respectively.
if the input was broadbanded, $X[k] \ne 0$ for all $k$.  so you can divide $Y[k]$ with $X[k]$.
the result: $$ H[k]=\frac{Y[k]}{X[k]} $$ is your system response.  this is called the "Dual-channel FFT" method of measuring system response. if there is a delay, you will see a linear trend to the unwrapped angle of $H[k]$. the negative of the slope of that linear trend is proportional to the delay.

you can also use cross-correlation between the input and output to discover the prominent delay of the system.  the definition of that is

$$ R_{xy}[m] = \sum_n x[n] y[n+m] $$

i'm deliberately playing fast-and-loose with the limits of the summation.  the value of $m$ that maximizes $R_{xy}[m]$ is a good candidate for the delay of the system.


Deriving difference equation I have this block-diagram of Karplus-Strong basic algorithm.

Where $$H_{d}(z) = \frac{1+z^{-1}}{2}$$
The problem occurs when I try to get the difference equation. As I can see from the diagram, that would be
$$y(n) = x(n-N)+\frac{y(n)+y(n-1)}{2}$$
but I've found online this one(which I know how to implement in code, unlike the first one)$$y(n) = x(n)+\frac{y(n-N)+y(n-(N+1))}{2}$$
Are these two equation somehow equivalent?

Neither of the two equations match with the diagram. It's actually quite straightforward: at the input of the delay you have

$$x[n]+\frac12\big(y[n]+y[n-1]\big)\tag{1}$$

Then a delay of $N$ samples is applied to compute the output, hence

$$y[n]=x[n-N]+\frac12\big(y[n-N]+y[n-1-N]\big)\tag{2}$$

But here I've found a different diagram, which corresponds to the last equation in your question.


recognizing frames as voiced/unvoiced/silence i have to implement a function in matlab,for labeling silenced,voiced and unvoiced frames in signal.
I mean i want to specify each frame as a voiced or unvoiced or silenced frame by using
pitch frequency and zero processing rate and energy.
how can i calculate pitch frequency?
where should i start from?
Thanks

I think that when you said "zero processing rate" you meant  "zero crossing rate", are a lot of way to find the pitch frequency, autocorrelation based methods is one choice, if your signal is always monophonic you can try AMDF (Average Magnitude Difference Function), AMDF is calculated  by the difference between the waveform summing a lagged version of itself, I like to use AMDF, work fine in monophonic signals, this can be defined as:

$$AMDF(m)=\sum_{i=0}^{n-1-m}|x(i) - x(i+m)|$$

Where:

x=your signal
n=length(x)-1
i = 0 ... n
m = time lag block size

Together Pitch Frequency + Zero Crossing Rate + Energy, can track if your frame is voiced/unvoiced/silenced

The pitch Frequency say if exist periodicity in your frame
Zero Crossing Rate say the times that amplitude of the signal changes
(crosses the x-axis)
Energy says how high or low is your frame

now you need to make decisions based in some threshold, and based in this decisions you can label your frames