Converting transfer function that is a sum of unusual rational polynomials to finite difference equation I have the following rather exotic transfer function:

$$
H(z) = cz^{-m} + \frac{b_0 z^{-1} + b_1z^{-2} + \dots + b_{2m}z^{-2m}}{1 + a z^{-1}} + \frac{q_0 z^{-1} + q_1z^{-2} + \dots + q_{2m}z^{-2m}}{1 + p z^{-1}}
$$

Where $c$, $b_0,\ldots,b_{2m}$,  $a$, $q_0,\ldots, q_{2m}$, $p$ are real valued constants and $m > 1$.

How do I go about converting this to a finite difference equation via the inverse $\mathcal Z$-transform?

I thought I had got to grips with the inverse $\mathcal Z$-transform, but I can't seem to convert $H(z)$ to a finite difference equation that produces the right output.

My attempt:

$$\begin{align}
y[n] &= c x[n - m] &+& (b_0 x[n - 1] + \dots + b_{2m} x[n - 2m]) \\
&&-& a y[n - 1]\\
&&+& (q_0 x[n - 1] + \dots + q_{2m} x[n - 2m])\\
&&-& p y[n - 1]
\end{align}$$

This doesn't seem to be correct. When I test this its magnitude response is not even close to the original transfer function.

My assumption is that the three terms $cz^{-m}$, $\frac{b_0 z^{-1} + b_1z^{-2} + \ldots + b_{2m}z^{-2m}}{1 + a z^{-1}}$, $\frac{q_0 z^{-1} + q_1z^{-2} + \ldots + q_{2m}z^{-2m}}{1 + p z^{-1}}$ can be transformed separately.


What is a filter bank? Suppose, I have 3 kernels:

$$ \left[
\begin{array}{cc}
  a&b&c\\
  d&e&f\\
  g&h&i
\end{array}
\right] $$
$$ \left[
\begin{array}{cc}
  p&q&r\\
  s&t&u\\
  v&w&x
\end{array}
\right] $$
$$ \left[
\begin{array}{cc}
  \alpha&\beta&\gamma\\
  \delta&\epsilon&\zeta\\
  \eta&\theta&\iota
\end{array}
\right] $$

How can I create a filter bank from them?

Should I AND, or OR, or add them together?

Or, should I just apply each of them one by one to my test-image through three separate convolution operations?

Since the term linear does not appear in the question and the current answers, let me offer a complementary perspective.

A kernel in this acceptation (especially for images, which don't always follow linear rules, think about occlusion or saturation) is an array that is applied, somehow, on any input data. One often distinguishes linear and non-linear kernel (as one has linear and nonlinear filters, even is the terminology might seem improper).

Let us start from the linear point-of-view in the most specific sense: the filter array is applied as a convolution. Then @MarcusMuller's answer is perfect: a set, an array of linear filters, applied to input data as  convolutions to yield several separate output data.
Most  additional scalar linear operation (like the sum, the average, a weighted combination) on the output would be "useless": as they commute, summing the output is equivalent to summing the three filters in one single filter, and performing only one single convolution on the data.

Which leads us back to the objective in your comment; traditionally, a linear (analysis, I will come back on that later) filter-bank (FB) is used to split or separate data into components, often with separate spectra or a narrower ocontent (low-, mid- or high-frequencies for a three-band filter-bank). Or to merge different data streams into  others, with a wider spectrum. So a generic multi-input-multi-output (MIMO) FB takes one or several inputs, filters them into one or several outputs. One then distinguishes analysis or synthesis filter banks.

Generally, recombining the outputs from an analysis FB drifts away from the separation objective. But a single filter is also a filter-bank (not very interesting per se though).
But sometimes, this can be more efficient (computationally for instance).

Now, having narrower/wider outputs invites to rate variations, like downsampling and upsampling before or after the filters. To me, the most accepted sense of a filter bank is a bank of linear filters optionally combined with (linear, but not shift-invariant) upsampling or downsampling operations. And it is somewhat related to linear transforms, allowing expansion or shrinkage on the number of coefficients (they can be critical, oversampled or undersampled).

Then, people extend the notion to nonlinearity: filters can be nonlinear (like the median) and the kernels are interpreted as weights applied to a chunk of data. Or the data can be combined in non-linear ways, with $\min$, $\max$, AND or OR...

But in your case, as Marcus said, I'd bet on three standard filtered outputs. But in this case, there is not relationship between filters (except their kernel size), and what is powerful in the filter-bank theory is the linkage between the filters, and how one can optimize them. Now a couple of pointers:

Filter bank (wikipedia)
What is a Filter Bank?
Multirate filter bank and multidimensional directional filter banks


Are the Daubechies 4 wavelet function and scaling function dimensionless? If the dimension of the independent variable $x$ of $f(x)$ is length, i.e. $[x]=L$, then what is the dimension of Daubechies 4 (hereafter D4) wavelet and scaling function?

Generally yes, discrete wavelets are often implemented as dimensionless filters, like discrete averagers or differences.
Additionally, wavelets are non dependent on the variable. If this is pertinent to your application, it is possible to equip the wavelet with dimensions: morally, the scaling coefficients ought to be of the same dimension as the variable, yet the detail coefficients can be thought as homogeneous to derivatives (at level 1), second derivatives (at level 2) etc. Consequently, with the proper normalization, they could be expressed in $\textrm{m}/\textrm{s}$, $\textrm{m}/\textrm{s}^2$, etc. according to each level.
If I detail the computation for the analysis/synthesis, decomposed onto the primal wavelet space ($\phi$ and $\psi$) and reconstructed in the dual wavelet space ($\tilde{\phi}$ and $\tilde{\psi}$):
$$f(x)=\sum <f,\phi_{j,k}>\tilde{\phi}_{j,k}(x) + \sum <f,\psi_{n,k}>\tilde{\psi}_{n,k}(x)\,.$$
The dimensional analysis entails that products or convolutions of $\phi$ and $\tilde{\phi}$, and of $\psi$ and $\tilde{\psi}$ should be unit-less. In other words, I would choose them  with "inverse" dimensions, and this could be done per scale.


Filtering Angular Measurements I am trying to apply a filter to reduce noise in angular measurements. Since the data comes with noise, the filter needs to be able to handle measurements that are not smoothly progressing; for example, the filter should be able to work with a series like below:
2, 358, 3, 5, 10, 18 ... 304, 315, 330, 348, 9, 15

As you can see, the second measurement of 358 is an erroneous measurement. How can I reliably put this into a linear scale so that I can apply a regular filter like a moving average filter to the data.
I was initially thinking of setting a maximum change parameter, like 90 degrees, and if such a change is observed to correct that.
I am very new to signal processing. Maybe this is a very common problem that is easily handled. Hope you  can put me in the right direction.

Median filtering is non-linear, and pretty awesome about removing outliers.  You just need to adjust the length of the filter based on the estimate of the frequency of the errored samples.


When inverting a transfer function, solving for the input using the output does the causality status change suppose $y(n)=ax(n-1)+bx(n-2)+\dots$  ($y$ is the output and $x$ the input). What happens if I want to solve $x(n)$ from $y(n)$?

Z transform: $$Y(z)=G(z)X(z)\tag{1}$$
then         $$X(z)=\frac{1}{G(z)}Y(z)\tag{2}$$

What are the properties of $1/G(z)$ ? If $(1)$ is causal what is the status of the inverse $(2)$? The roles  of the poles and zeros have changed.

As you've pointed out, inversion leads to poles at locations of the zeros of the original transfer function and vice versa. Assuming that $G(z)$ is causal and stable (i.e., it has all its poles inside the unit circle), we have to distinguish $3$ cases:

$G(z)$ has at least some zeros outside the unit circle. This means its inverse has some poles outside the unit circle, and consequently, it cannot be causal and stable. If $G(z)$ has no zeros on the unit circle, there exists a stable impulse response corresponding to $1/G(z)$ but it cannot be causal. This is because a transfer function does not uniquely determine an impulse response. We can get different impulse responses corresponding to different regions of convergence of $1/G(z)$.
$G(z)$ has some of its zeros on the unit circle. No stable inverse exists because $1/G(z)$ has poles on the unit circle.
$G(z)$ has all its zeros inside the unit circle, i.e., it is a minimum-phase system. Consequently, $1/G(z)$ also has all its poles and zeros inside the unit circle (i.e., it is minimum-phase), and can be implemented by a causal and stable system.


How to determine Hilbert coefficients I want to take a real signal on my FPGA and create a complex one.  I have FIR core implemted for a Hilbert Transform, and Xilinx's datasheet on it shows coefficients of: (-819,0,-1365,0,-4096,0,4096,0,1365,0,819) as an example.

I wanted to add more taps, so as practice I used the matlab fdesign function:

d=fdesign.hilbert(11,0.1);
hd=design(d,'firls');
coef=hd.Numerator

The coefficients seem reasonable (though I will need to multiply by 10000 to get integers).  What I don't understand is the 0.1 transition width variable (I took the 0.1 from an example online), nor if I need to tune this filter for a frequency range, or if they are generic enough.

The resulting analytic signal after using a the Hilbert transform filter is ideally a signal without any negative frequencies, hence, attenuating -1 (using normalized frequency) to 0 while transmitting 0 to 1 (note that this is a complex output and therefore one do not have the standard symmetry around 0). This is obtained by phase-shifting positive frequencies pi/2 and negative frequencies -pi/2 while keeping the magnitude to be one and adding the output with the original signal.

Since it is not possible to obtain this change in phase-shift without a transition band one will need to specify how wide the transition band should be.. Hence, in your case the 0.1 would indicate that your designed filter will perform well from 0.05 to 0.95 with a transition band from -0.05 to 0.05 and similarly from 0.95 to 1.05.

Increasing the order with constant transition band will produce a higher attenuation/smaller approximation error in the passband, while decreasing the transition band with the same order will have the opposite effect.

Also, if you will implement it on an FPGA I would suggest multiplying the coefficients with a power of two, e.g. 8192 to get a better range for the output and use of the bits.


Why are the mark and space tones of these FSK signals reversed? I'm working on a project involving FSK-modulated signals. I'm trying to reverse-engineer these signals to determine the preamble pattern used, because the protocol is proprietary.

Shown below is a plot of the real parts of two captured transmissions that seem to follow the same pattern - except their frequencies appear to be inverted. This is further evidenced by the instantaneous frequency plots also shown. I captured these signals using an RTL-SDR dongle, which outputs complex I/Q samples. What does this mean?

The apparent exchange of preamble symbols between mark and space in separate transmissions indicates that the packet is differentially encoded in some way.  NRZI encoding is used in both Packet Radio and AIS, so you might wish to investigate if that appears to be the case for your unknown packet format as well.


Time and Frequency Domain Relationship Why do we say that a good time response means a bad frequency response?

These terms refer to the general duality properties of these two domains. It should be obvious that higher frequency components mean faster changes in time and lower rise/fall times. In most natural systems, some of these fundamental properties oppose each other:
Shorter rise time means more overshoot and higher settling times and vice versa; lower overshoot means a higher settling time as well as a higher damping factor resulting in a higher rise time. In general, achieving better frequency domain properties (in the sense of time domain duals) means worsening other time domain properties. Apart from that, in specific applications such as digital transmission, cleaner pulses (lower rise/fall time) means broader spectrum and frequency content, which is undesirable in the sense of bandwidth efficiency. But again there is no unique answer to what you have asked. Each application must consider its own tradeoffs.


When should I use a frequency spectrum rather than a spectrogram representation? It seems to me that the spectrogram representation of a finite signal gives you significantly more useful information than a frequency spectrum.  As I understand it, the spectrogram breaks down the spectrum over time whereas the frequency spectrum gives you no time dependent information.  Why wouldn't I want to use this kind of representation all the time?  When would it be more useful to use the frequency spectrum rather than the spectrogram?  Is it just computationally expedient to produce a frequency spectrum?

The spectrogram suffers from a time/resolution restriction - the smaller your time step, the smaller your analysis windows become, the worse your frequency resolution gets.

It is a waste of frequency resolution to break down a stationary signal into tiny temporal slices. Thus, if you know in advance that your signal you analyze is stationary; or if you have pre-segmented your signal into stationary segments (using a transient/onset detector; novelty detection; or an auxiliary signal/sensor), there is no reason to lose frequency resolution by using the spectrogram.


Find Fourier Transform of signal that satisfies this equation I have attempted it this way:

$G(j\omega)$ is a rectangular pulse in frequency domain. Using results, inverse Fourier transform of $G(\omega)$ will be $$g(t)=\frac{\sin(2t)}{\pi t}$$

Now given that $g(t)=x(t)\cos(t)$

Using $\sin (2t)=2\cos(t)\sin(t)$, I got $$x(t)=2\frac{\sin(t)}{\pi t}$$

But, I am struggling to draw the parallel for part b of the question.

I was using the multiplication property for part b, and it turned out,

$$G(j\omega)=\frac{X\left(\omega-\frac 23\right)+X\left(\omega+\frac 23\right)}{2}$$

How to find $X(\omega)$ and proceed from here onwards?

Your procedure for the first exercise is correct. Note however that, as you may have noticed already, that approach is only useful if the cosine has unitary pulsatance. There is a generic way to do this, if we think of the problem in the frequency domain.

I'll apply this method to the exercise you already solved, so that you can work the second one on your own. In the frequency domain, your exercise can be stated as follows:

$$G(\omega)=\frac{1}{2\pi}X(\omega) *[\pi\delta(\omega-1)+\pi\delta(\omega+1)]$$

Remember that convolving signals with Dirac deltas can be easily made graphically. It's just centering the signal being convolved (in this case, $X(\omega$)) in the position where the impulses are.

Notice that when the impulses are centered in $\omega=\pm1$, if $X(\omega)$ was a rectangular pulse of width $2$ and amplitude $2$ as well, then the convolution would return a rectangular pulse of width $4$ (one for each impulse, and as they don't overlap they keep the rectangular shape) of amplitude $2\pi$. But because convolution in the frequency domain has to be divided by $2\pi$ when going back to the time domain, we get that if $X(\omega)=\mathrm{rect}(\omega/2)$, then the equality holds. That transform corresponds indeed to the signal $x(t)$ you found using the trigonometric identity.

See that this approach is useful for any position of the Dirac deltas. In the second case, they are at $\omega=\pm 2/3$. Can you solve it now?


Compressive Sensing - Incoherence Property Compressive Sensing is built on 2 properties: 1) the sparsity of the representation basis relative to the sampling basis and 2) the incoherence between the singular vectors from each of the 2 bases in a). On the surface this seems fine to me but he "incoherence" relationship is confusing me a little.

Some texts refer to the coherence between the bases (representation, sampling) and other refer to the coherence between each basis and the standard basis (e1, e2, e3 ...).

Is there a difference between these two statements?

  "We are in the position to state our main result: if a matrix has row
  and column spaces that are incoherent with the standard basis, then
  nuclear norm minimization can recover this matrix from a random
  sampling of a small number of entries."

Page 6: http://statweb.stanford.edu/~candes/papers/MatrixCompletion.pdf

and

"Incoherent sampling ... The coherence between the sensing basis and and the representation basis".
Page 3: http://authors.library.caltech.edu/10092/1/CANieeespm08.pdf

My question is related to this question:
https://dsp.stackexchange.com/a/13017/4038

I am trying to answer your question about incoherence here rather than update my previous answer on another question of yours.

Compressive sensing requires low coherent pairs. So the lower $\mu(\Phi,\Psi)$, the better. Actually is $\Phi$ is spike basis (identity matrix) with $\phi_k(t) = \delta(t-k)$, and $\Psi$ is Fourier basis with $\psi_j(t) = 1/\sqrt n e^{-i \cdot 2\pi \cdot jt/n}$, which is just the example I showed you in another question, $\mu(\Phi,\Psi) = 1$, and the maximal incoherence is achieved. It is similar that if you apply the spike basis and the orthogonal basis obtained from SVD.

Other coherence includes but not limited to: coherence between noiselets and Haar wavelets is $\sqrt2$; coherence between noiselets and Baubechies D4 and D8 are $2.2$ and $2.9$, respectively; random matrices are largely incoherent with any fixed basis $\Psi$.

You may refer to Candes-Romberg theorem for more details on incoherent sampling.


Matlab stores all level of Wavelet's Approximation & Details in a single Row I'd like to ask when we apply Wavelet transform on 1-D Signal in Matlab it simply stores both Approximations and Details into a single row which is quite annoying for me,

have a look,

Where,

X = original Signal
cDn = Details coefficients
cAn = Approximations

This is how Matlab do it,

Suppose I have a random sequence,

X=[2 2 3 100 4 0 98 100 90 2 3 67 98 0 6 6 89 9 21 78]

If I apply 2-level Wavelet transform on this signal it will be decomposed in to cA1 and cD1. cA1 will then further decomposed into cA2 and cD2. But matlab will store all these A's and D's into a single row. Have a look what matlab will do to the resultant coefficients after performing Wavelet Transform on it.

X = 20 Samplesc  -> cD1 = 10 Samples
                    cA1 = 10 Samples
                                     -> cD2 = 5 Samples
                                        cA2 = 5 Samples

Rows = 1
Columns = 20

Column 1 to 5 =   cA2
Column 6 to 10 =  cD2
Column 11 to 20 = cD1

Why Matlab does that ? Why it doesn't decompose signal in Rows since that would be easier to represent.

You can use the function dwt instead that does single-level decomposition and use that again for 2-level decomposition. Personally, I use wavedec all the time, I actually find it easy to have it in one array. Boils down to personal preferences I guess.


Minimum FFT sample time - guitar frequencies and harmonics I am working on a project whereby I am trying to determine the exact note played on an electric guitar with a microcontroller (as quickly as possible after the string is plucked) and to use this information (ie digital value of note) to manipulate an audio output via DAC.

I am using an STM32F4, sampling the DC offset and active fitered guitar signal with a 12-bit ADC (sample rates pretty much as fast as I need with timer triggered ADC DMA) and the highest number of F32 FFT points I can achieve with the RAM available is 4096 (ie 8192 samples).

The guitar frequencies I need to measure are between 82Hz and about 831Hz, with a resolution at the low end of 2Hz and about 22Hz at the top.

Am I correct in saying minimum sample time to achieve this is about 0.5s, driven by the minimum frequency and resolution requirement?
And at higher frequencies I can sample for shorter (with a higher sample rate)?
Or is there (hopefully) something I'm missing?
Can I measure the lower frequencies with fewer FFT points (and thus samples), saving time?

Your highest frequency is 831 Hz; since undersampling to ignore the sub-82 Hz would only save you 10% sampling rate, let's ignore that possibility.

So, Nyquist says your sampling rate should be at least 1.662 kHz; now, your problem is, as @fibonatic mentioned, that a the freq. resolution of any DFT is alwas

$$\Delta f = \frac{f_\text{Nyquist}}{N_\text{DFT}}\\
\implies\\
N_\text{DFT} = \frac{f_\text{Nyquist}}{\Delta f}$$

and the time taken by capturing of samples you'll need for that DFT is

$$T = \frac{N_\text{DFT}}{f_\text{sample}} = \frac{N_\text{DFT}}{2f_\text{Nyquist}} = \frac{\frac{f_\text{Nyquist}}{\Delta f}}{2f_\text{Nyquist}}= \frac{1}{2\Delta f}$$

So no matter at which rate you'd sample, there's nothing oyu can do about the time period needed for a full-resolution DFT.

You could just pad a shorter observation with zeros, and that will increase your resolution. Point is that you trade resolution for quality, and the variance is going to grow significantly.

The DFT is simply not the right choice here; if you need 2Hz resolution at the lower frequency, you get the same at the higher frequencies without any need.

Honestly, think about it: you're planning to use a DFT that gives you 4096 frequency bins, where you'd need to distinguish maybe 100 tones, if that many.

I'm pretty certain that the cheap guitar tuners you get don't work that way. I'd assume they simply employ dedicated hardware (not firmware like you'd need to write) to implement a simple filterbank – possibly even a IIR filterbank, and if I had a guess, at very low sample bit width, with simple biquads. Or, even simpler, they just have couple of analog filter bandpasses, and just undersample the resulting sample streams and do the same digital filters on all the octaves they preselect that way.

If I had to do this, I'd atually do the same on your microcontroller (by the way, 8192-FFT sounds gigantic for a Cortex-M4, but you know your device better than I do), but with different architecture. Do a rough calculation: If you keep the sampling rate low, how many FIR filters of the different lengths you'll need can you let run over the data? Does implementing your filterbank as a PFB filterbank help? Does making a hierarchic decision (e.g. just first an "octave selector", then depending on the result of that, throwing away half, $\frac38$,$\frac78$,… of the samples and then doing a tone-detection filterbank.


Cepstrum Calculation of Rational Function H(z) I am trying to solve my first problems at cepstrum calculation. I want to calculate the complex cepstrum $\hat{h}[n]$ of a signal $h[n]$ with Z-Transform:
$$H(z)=\frac{(1-0.5z^{-1})(1+4z^{-2})}{(1-0.64z^{-2})}$$
I know that $\hat{H}[z]=\log(H(z))$ and $\hat{h}[n]=\frac{1}{2π}\int_{-π}^{π}\log[H(e^{jω})]e^{jωn}dω$

Which is the easiest way to approach this?

I started by applying $\log()$ to H(z) which gave me:
$$\log(H(z))=\log(1-0.5z^{-1})+\log(1+4z^{-2})-\log(1-0.64z^{-2})\Rightarrow $$
$$\log(H(z))=\log(1-0.5z^{-1})+\log(4(z^{-2}+\frac{1}{4}))-\log((1-0.8z^{-1})(1+0.8z^{-1}))\Rightarrow$$
$$\log(H(z))=\log(1-0.5z^{-1})+\log(4(z^{-1}+\frac{j}{2})(z^{-1}-\frac{j}{2}))-\log((1-0.8z^{-1})(1+0.8z^{-1}))\Rightarrow$$
$$\log(H(z))=\log(1-0.5z^{-1})+\log(4)+\log(z^{-1}+\frac{j}{2})+\log(z^{-1}-\frac{j}{2}))-\log(1-0.8z^{-1})-\log(1+0.8z^{-1})$$

Is this apporach right? If so, how do I go on? Any help is appreciated. Thanks in advance!

We're looking for a series representation

$$\log H(z)=C(z)=\sum_{n=-\infty}^{\infty}c[n]z^{-n}\tag{1}$$

which converges in an annular region $r_1<|z|<r_2$ with $0<r_1<1<r_2$.

Note that the poles as well as the zeros of $H(z)$ lead to singularities of $C(z)$. All poles and zeros of $H(z)$ inside the unit circle contribute to the right-sided part of $c[n]$, whereas the poles and zeros of $H(z)$ outside the unit circle contribute to the left-sided part of $c[n]$. Consequently, a minimum-phase transfer function $H(z)$ has a right-sided cepstrum.

The given function $H(z)$ has all its poles inside the unit circle ($z^{\infty}_{1,2}=\pm 0.8$, $z^{\infty}_{3}=0$), one of its zeros is inside the unit circle ($z^{0}_{1}=0.5$), and two zeros are outside the unit circle ($z^{0}_{2,3}=\pm 2j$). Consequently, its complex cepstrum $c[n]$ is a two-sided sequence.

In order to find the sequence $c[n]$ we make use of the well-known Mercator series:

$$\log(1-z)=-\sum_{n=1}^{\infty}\frac{z^n}{n},\qquad |z|<1\tag{2}$$

From $(2)$ it is clear that terms of the form

$$\log(1-az^{-1})=-\sum_{n=1}^{\infty}\frac{a^n}{n}z^{-n},\qquad |z|>|a|\tag{3}$$

contribute to the right-sided part of $c[n]$, whereas terms of the form

$$\log(1-bz)=-\sum_{n=1}^{\infty}\frac{b^n}{n}z^{n}=\sum_{n=-\infty}^{-1}\frac{b^{-n}}{n}z^{-n},\qquad |z|<\frac{1}{|b|}\tag{4}$$

contribute to the left-sided part of $c[n]$. Note the argument can be easily generalized to higher powers of $z$ on the left-hand sides of $(3)$ and $(4)$.

The term $(1+4z^{-2})$ in the numerator of the given function $H(z)$ corresponds to two zeros outside the unit circle, and, consequently, it contributes to the left-sided part of $c[n]$. Hence, according to $(4)$, it must be rewritten in terms of positive powers of $z$:

$$(1+4z^{-2})=4z^{-2}(1+0.25z^2)\tag{5}$$

We have

$$H(z)=z^{-2}\frac{4(1-0.5z^{-1})(1+0.25z^2)}{(1-0.64z^{-2})}=z^{-2}\tilde{H}(z)\tag{6}$$

I think that due to the linear-phase term $z^{-2}$ in $(6)$ we actually cannot compute the complex cepstrum of $H(z)$, but we can compute the complex cepstrum of $\tilde{H}(z)$, which is just a shifted version of the original function $H(z)$.

From $(6)$ we obtain

$$\log \tilde{H}(z)=\log(4)+\log(1-0.5z^{-1})+\log(1+0.25z^{2})-\log(1-0.64z^{-2})\tag{7}$$

Using $(3)$ and $(4)$ we get

$$\log \tilde{H}(z)=\log(4)-\sum_{n=1}^{\infty}\frac{(0.5)^n}{n}z^{-n}+\sum_{n=-\infty}^{-1}\frac{(-0.25)^{-n}}{n}z^{-2n}+\sum_{n=1}^{\infty}\frac{(0.64)^n}{n}z^{-2n}\tag{8}$$

which is equivalent to

$$\log \tilde{H}(z)=\log(4)-\sum_{n=1}^{\infty}\frac{(0.5)^n}{n}z^{-n}+\sum_{n=-\infty\\n\textrm{ even}}^{-1}\frac{2(-0.5)^{-n}}{n}z^{-n}+\sum_{n=1\\n\textrm{ even}}^{\infty}\frac{2(0.8)^n}{n}z^{-n}\tag{9}$$

The region of convergence of the series representation $(9)$ is $0.8<|z|<2$, which includes the unit circle, as required by the definition of the complex cepstrum.

From $(9)$ we can directly write down the complex cepstrum of $\tilde{H}(z)$:

$$\tilde{c}[n]=\begin{cases}\log(4),&n=0\\
\frac{1}{n}\left[2(0.8)^n-(0.5)^n\right],&n>0 \land n\textrm{ even}\\-\frac{1}{n}(0.5)^n,&n>0 \land n\textrm{ odd}\\\frac{2}{n}(-0.5)^{-n},&n<0 \land n\textrm{ even}\\0,&n<0\land n \textrm{ odd}\end{cases}\tag{10}$$


Complex time domain signals and QAM Is there an inherent connection between complex time domain signals and QAM?

Given a complex time-domain signal is it always implied that this is a QAM-demodulated form of a real time domain signal?

Are there any other commonly used ways of converting a real time domain signal to a complex time domain signal that don't imply QAM?

Given a complex time-domain signal is it always implied that this is a QAM-demodulated form of a real time domain signal?

No, we generally use complex baseband to describe all kinds of signals, or even passband systems (e.g. frequency-selective or time-variant channels). So, any bandlimited signal has an equivalent complex baseband, and it's not implied this has anything to do with QAM or even digital communication.

QAM is just one example of something that is really easy to write down in complex baseband.

  Are there any other commonly used ways of converting a real time domain signal to a complex time domain signal that don't imply QAM?

I don't even know how QAM comes into this – you transform any real-valued bandpass signal to complex baseband through multiplication with its (complex) carrier and low-pass filtering. QAM really has nothing to do with it.


Discretize process noise in Kalman filter Reading P. Andrews et al. I see that it is very common to do the following approximation of the process noise covariance matrix:

$$Q_{k} = G_{k-1}QG_{k-1}^{T}\Delta t$$

so that the propagation becomes:

$$P_{k} = F_{k-1}P_{k-1}F_{k-1}^{T} + Q_{k}$$

however, when I try it in a simulation I need to scale it once more:

$$P_{k} = F_{k-1}P_{k-1}F_{k-1}^{T} + Q_{k}\Delta t$$

else the discretization is wrong - I am not sure why this happens.

The simulated process noise is by adding the following to the ground truth:

velocity_measured = velocity_true + sigma * np.random.randn(1)
position = position + velocity_measured * dt

so that the process noise before discretization becomes:

Q = sigma * sigma.

Hi: I've never dealt with discretization in the kalman filter ( my models were already discrete ) so take the following  with a level of uncertainty ( no pun intended ). Also, you didn't show the original equations so I'll refer to them as the observation equation and the system equation. Based on your updating equations,  it seems that you're using the non-lagged formulation where the observation equation and the system equations are updated at the "same" time rather than the observation equation being one lag ahead.

So, assuming the above is correct ( without the original equations, it's a little hard to tell ),  the reason it's multiplied by $\triangle t$ twice is because there are two updates to the system noise due to discretization of the two equation process.  The first update is due to the $G_k$ matrix multiplying the variables in the observation equation at time $k$.  That's what causes the first $\triangle t$ multiplication. But then, you also have the update of the system equation at time $k$ due to the multiplication of the state variables by $F_{k-1}$. This is why
$\triangle t$ is used a second time.

But, rather than think of the discretization as using multiplcation twice, you can think of it as using multiplication sequentially. It looks as though once is multiplying twice because, in your formulation, the observation equation and system equation are not lagged with respect to each other. But it's really just sequential multiplication due to sequential updating. Note that if you used the other KF formulation where there is lagging, then the discretization process would probably be easier to see. I hope this helps some.

                                                 Mark


Discrete-time system: getting the true frequency from normalized frequency Lеt's say that i have a plot where i can see the spectrum of some discrete signal and let's say that it's x axis is normalized (by normalized i mean $\omega/\pi$ on the $x$ axis). Now, i know that in order to get to the real frequency at some particular point on the plot, i need to know the sampling frequency $f_s$ and then multiply that particular point with $f_s/2$ in order to get the actual frequency.

Now, i know this is maybe something trivial, but i am having struggle to understand why is this so, i mean, why multiplication by $f_s/2$ gives us the actual frequency in hertz?

Any help appreciated!

The maximum frequency of a discrete-time signal is half the sampling frequency. A signal with maximum frequency in discrete time is an alternating sequence, and since its period is two sample intervals, its frequency is $f_s/2$. All frequencies in the discrete domain are normalized by that maximum frequency.

This implies that a discrete-time system's behavior can be described independently of the sampling frequency. E.g., if you have a low pass filter, the cut-off frequency is fixed in terms of relative frequency. The actual cut-off frequency in Hertz depends on the sampling frequency.

Often you will encounter the normalized frequency in radians:

$$\omega=2\pi\frac{f}{f_s}\tag{1}$$

where $f_s$ is the sampling frequency, and $f$ is the frequency in Hertz. Note that the frequency response of a discrete-time system is $2\pi$-periodic with respect to $\omega$, i.e., with respect to $f$ it is periodic with period $f_s$.


Determine Phase Difference and Time Delay Between Two Pulsed Sinusoidal Signals I am working on a code that generates two signals: A pulsed sinusoidal signal ($1\textrm{ MHz}$) beginning at $0$ phase, and lasts for say $5\ \mu\text{s}$, then trails zeros. The other signal has the same pulse width and frequency ($5\ \mu\text{s}, 1\textrm{ MHz}$), but begins at a different time ($1\ \mu\text{s}$ later for example), and with a different phase ( $\pi$ radians).

If either the phase or time delay were consistent between the two, then I would be able to solve this problem, but how do you go about solving this when both variables differ between the signals?  Any feedback on this is greatly appreciated.

Pulse width and frequency alone are insufficient information to determine phase (relative to a fixed point in time or to each other) and delay.

The main difference between two equal length sinusoids with different starting phases is in the difference between their starting and ending transients.  So after estimating the time difference between the two pulse envelopes, you would have to look at how each pulse envelope starts and ends (which requires having a lot more bandwidth than narrow band frequency estimation alone needs).


Python tool for time-frequency analysis I am trying to perform time-frequency analyses using the PyWavelets (pywt) toolkit for python. My ultimate goal is to perform time-frequency analyses for EEG signals but I am starting with something simpler.
For a sanity test, I am creating a simple signal of length 2 seconds, with sample rate 250Hz, containing 2 sine waves - one of 3Hz and one of 10Hz. I would like to create a time-frequency plot that has two horizontal lines - one for the 3Hz and one for the 10Hz, which looks like this (only for illustration purposes):

For this purpose, I tried using code from the following tutorial : http://ataspinar.com/2018/12/21/a-guide-for-using-the-wavelet-transform-in-machine-learning/, specifically in section 3.1 of the tutorial.

This is a minimal example based on the code from the tutorial:

from UliEngineering.SignalProcessing.Simulation import sine_wave
import pywt
import numpy as np
import matplotlib.pyplot as plt

def plot_wavelet(time, signal, scales,
             waveletname='cmor',
             cmap=plt.cm.seismic,
             title='Wavelet Transform (Power Spectrum) of signal',
             ylabel='Period (seconds)',
             xlabel='Time'):
    dt = time[1] - time[0]
    [coefficients, frequencies] = pywt.cwt(signal, scales, waveletname, dt)
    power = (abs(coefficients)) ** 2
    period = 1. / frequencies
    levels = [0.0625, 0.125, 0.25, 0.5, 1, 2, 4, 8]
    contourlevels = np.log2(levels)

    fig, ax = plt.subplots(figsize=(15, 10))
    im = ax.contourf(time, np.log2(period), np.log2(power), contourlevels,
    extend='both', cmap=cmap)

    ax.set_title(title, fontsize=20)
    ax.set_ylabel(ylabel, fontsize=18)
    ax.set_xlabel(xlabel, fontsize=18)

    yticks = 2 ** np.arange(np.ceil(np.log2(period.min())),
    np.ceil(np.log2(period.max())))
    ax.set_yticks(np.log2(yticks))
    ax.set_yticklabels(yticks)
    ax.invert_yaxis()
    ylim = ax.get_ylim()
    ax.set_ylim(ylim[0], -1)

    cbar_ax = fig.add_axes([0.95, 0.5, 0.03, 0.25])
    fig.colorbar(im, cax=cbar_ax, orientation="vertical")
    plt.show()

def generate_sine_wave(length, samplerate, frequencies):
    wave = np.zeros(int(length * samplerate))
    for frequency in frequencies:
        wave += sine_wave(frequency=frequency, samplerate=samplerate,
        length=length)
    return wave

signal = generate_sine_wave(2, 250, [3, 10])
N = len(signal)
t0 = 0
dt = 1/250
time = np.arange(0, N) * dt +t0

scales = np.arange(1, 256)
plot_wavelet(time, signal, scales)

This plot from this code doesn't give me the plot I want, it looks like this:

I tried many modifications for this code but none gave me the result I want. And there are a couple of things I don't understand about the code:
- What is the purpose of the "period" variable in the "plot_wavelet" function and how do I make the y-axis show frequencies instead?
- What is the purpose of the "scales" variable?
- How do I define a frequency range that I want the result to include?
- How do I use linear scaling for the frequencies instead of log scale?

If anyone can give some pointers regarding this I will be very happy. Been spending some time trying to plot normal time-frequency plots but still haven't been able to find a python tool that performs this simple plot which makes sense to me.
Thank you,
Elad

You can find a nice tutorial for time-frequency analysis in Numerical python by Johansson, chapter 17. link to github repository.

You can also check the scipy.signal.spectrogram.

import numpy as np
from scipy import signal
from scipy.fft import fftshift
import matplotlib.pyplot as plt

# Generate a test signal, a 2 Vrms sine wave whose frequency
# is slowly modulated around 3kHz, corrupted by white noise
# of exponentially decreasing magnitude sampled at 10 kHz.

fs = 1e4
N = 1e5
amp = 2 * np.sqrt(2)
noise_power = 0.01 * fs / 2
time = np.arange(N) / float(fs)
mod = 500 * np.cos(2 * np.pi * 0.25 * time)
carrier = amp * np.sin(2 * np.pi * 3e3 * time + mod)
noise = np.random.normal(scale=np.sqrt(noise_power), size=time.shape)
noise *= np.exp(-time / 5)
x = carrier + noise

fig, ax = plt.subplots(2, figsize=(8, 7))

f, t, Sxx = signal.spectrogram(x, fs)
ax[0].pcolormesh(t, f, Sxx)
ax[0].set_ylabel('Frequency [Hz]')
ax[0].set_xlabel('Time [sec]')

# Note, if using output that is not one sided, then use the following:
f, t, Sxx = signal.spectrogram(x, fs, return_onesided=False)
ax[1].pcolormesh(t, fftshift(f), fftshift(Sxx, axes=0))
ax[1].set_ylabel('Frequency [Hz]')
ax[1].set_xlabel('Time [sec]')
plt.savefig("fig.png")
plt.show()


What is the standard way to handle images of different sizes for classification? I'd like to create an image classifier that takes a set of labeled images and creates a classifier to be used on some unknown images. Most of the examples I see assume all images are the same size (m x n pixels), but my training set will have pictures of various sizes. This is a problem because training my classifier will have the following steps:

Load each image into m x n array
Employ some filter
Transform m x n array into 1 x (m * n) array
Feature selection
Train classifier

The different sized images is a problem because step 3 will yield arrays of different lengths if the images don't have the same number of pixels.

I have a background in statistics, not image processing, and I haven't been able to find any general guidelines for how to handle this problem. I'll be working in Python, and it would also be nice to know if there are libraries perform image resizing.

Not sure about your specific application, but the Viola Jones face detection algorithm successively scans the candidate picture for faces in all possible sizes.

I.e., the algorithm search space is both the object (face) location and its size. The beauty of VJ's "cascade architecture" is that it rejects most hypotheses quickly, thus making the search feasible in real-time.


Find the $\mathcal Z$-transform of this function? I need to find the $\mathcal Z$-transform of $x(n)=a^{-n} u(n)$. Assume, $a$ is a positive constant , but the power of $a$ is negative.

Looking at the transform table, I found that $\mathcal Z$-transform for $x(n)=a^{-n} u(n)$ is not there, so I tried to found it by on my own, but I am not sure if my answer is right, please can anyone check it and correct me, in case?

Doing some calculations is appreciated. However, it is preferred that you write them in the question because a picture can vanish more easily. Do not hesitate to edit the math.

To help you with an hint in checking your answer, remember that:

$$ a^{-n} = \left(\frac{1}{a}\right)^n$$

a standard trick to play with known tables.


Calculate Impulse Response of a System From Its Input/Output Waveform via FFT From the theory, I knew something like :

$y(n) = x(n)*h(n)$

$Y(n) = X(n)H(n)$

$H(n) = \frac{Y(n)}{X(n)}$

Now I have $y(n)$ and $x(n)$ from a system, and $y(n)$ is 2176 samples long, $x(n)$ is 2176 samples long.

However, if I used following code to calculate impluse response, it wasn't correct to me, anyone could help figure out where I am wrong ?

bit = load('C:\Users\wangyang\Desktop\bit_wf.txt')

x = bit(:,2)

subplot(2,1,1)
plot(x)
title('x(n)')

sbr = load('C:\Users\wangyang\Desktop\wf.txt')

y = sbr(:,2)

subplot(2,1,2)
plot(y)
title('y(n)')

X = fft(x)
Y = fft(y)
h = ifft(Y/X)
plot(h) %plot impulse response

If I used ./ instead of / to calucalte the h as folloiwng:

X = fft(x)
Y = fft(y)
h = ifft(Y./X)
plot(h) %plot impulse response

The h looks like as following :

Thanks so much.

It would be best to upload the exact data for $x[n]$ and $y[n]$, but the reason I think of, is that either at least one of your signals is not a perfect baseband signal and has considerable amount of energy everywhere, or (this one is more probable) the signals are baseband but your sampling rate is not high enough, hence the Nyquist criterion is not satisfied and aliasing has occurred. In the latter case, try increasing the number of samples per each signal. Note that this is also an approximate, because the DFT package of MATLAB operates on discrete signals, yet with truncation in both time and frequency.


How to save Fourier series approximated signal to a WAV file I changed this Matlab/Octave code to approximate square wave by using combination of Fourier series and Fejér taper:
    % needed for Octave ------------------------------
    pkg load signal
    % ------------------------------------------------

    % Stage 1 ----------------------------------------
    figure
    N=20;                               % harmonics
    fs = 44100;
    fc = 440;
    w = 2 * pi * fc;
    r_start=0;
    r_end=10;
    r_step=1/fs;
    t = r_start:r_step:r_end;           % range
    w_sqr = square(w*t);          % Square wave [-1:1]

    plot(t,w_sqr);
    axis([r_start r_end -1.2 1.2]);
    grid on;
    hold on;
    audiowrite("squarewave.vaw", w_sqr, fs);

    % Approximation: Combined
    % Fourier series-Fejér taper (Tukey window?) -----------
    %i=1;
    %sum=0;
    %r=t;
    %for t=r
    %  for n=1:N
    %    %sum = sum + (2*sin(n*t)+sin(pi*n-n*t)-sin(n*t+pi*n))/(2*pi*n) % Fourier series
    %    %sum = sum + ((1-n/(N+1))*(-1)*sin(n*t)*(cos(n*pi)-1))/(n*pi); % Fejér taper
    %    sum = sum + ((n-2*N-2)*(cos(pi*n)-1)*sin(n*w*t))/(2*pi*n*(N + 1)); % (Fourier+Fejér)/2
    %  end
    %   Ffej(i)= sign(sum)*(1/2) + sum;
    %  i=i+1;
    %  sum=0;
    %end;
    %Ffej=Ffej';

    sum_ = sum((((1:N)-2*N-2).*(cos(pi*(1:N))-1).*sin((1:N).*w.*(0:1/fs:1)'))./(2*pi*(1:N)*(N + 1)), 2);
    Ffej = 2*sum_;
    Ffej = Ffej';

    plot(t,Ffej);
    axis([r_start r_end -1.2 1.2]);
    legend('Square', 'Fourier series-Fejér taper');

    audiowrite("approximated_squarewave.wav, Ffej, fs);

which results this plot:

I would like to save these signals to .wav files by using Matlab/Octave audiowrite() command.
What changes are needed to be done in code to save an approximation of 440Hz square wave sound sampled at 44100 Hz? My thought is that range t must be changed to support sample rate range and frequency but, how is the situation with calculation of Ffej ... or is data of Ffej directly usable in audiowrite() command?

OK, after digging the case by myself I managed to generate square wave files for both cases. Updated the source code found in my opening post.
In case of combined Fourier Series and Fejér taper approximations of square wave signal, result looks like this (plot done with Audacity):

and spectrum like this (Plot done with Octave):

(green = w_sqr, black dotted = Ffej, fc=1230Hz, fs=44100Hz, w_sqr and Ffej are data holders from source code)


Complex Morlet function and DWT I have an application that I must use a complex Morlet wavelet function (cmorfb-fc), I can run CWT (continuous wavelet transform) and it's fine, but I want to run DWT(Discrete wavelet transform) with this complex Morlet(cmorfb-fc) function too. then I want to compare the DWT and CWT results.
so Does anyone can help and explain, how can I implement DWT of complex Morlet (with the known Mallat filters)?

A continuous wavelet frame can be discretized with perfect reconstruction under some technical conditions. If this is possible, you obtain a discrete wavelet frame, often with some amount of  redundancy (depending on the wavelet shape). This is a result formalized by I. Daubechies, for instance.

For a long time, apart from the Haar or the Franklin wavelets, little was known about the possibility of obtaining a critical (non-redundant) discrete wavelet transform. Until the inception of the Meyer wavelet, and the subsequent ones (Daubechies, Symmlets, Coiflets, etc.). Their somewhat weird shapes is the sign that not every wavelet shape is possible. This is even truer as the complex Morlet is complex, hence somewhat already twice redundant for real signals.

So, Morlet cannot be critically discretized with 2-band (dyadic), as you already guessed, as you wanted to compare CWT and DWT. However, multi-band  filter banks are an option to build orthogonal counterparts of "Morlet"-like design. Malvar wavelets, Modulated complex lapped transforms or Lapped orthogonal transforms are  instances (see a Panorama of wavelets for references), and their can be  further decomposed across scales as $M$-band wavelets.

Or you can accept a non perfect reconstruction dyadic decomposition akin to Morlet.


Determining Causality and Time-Invariance of a system Consider the following system:
$$y(t-1)=\int_{-\infty}^\infty x()u(-t) d $$
where $u(t)$ is the unit step function, which is zero for $t<0$ and equals $1$ for $t>0$.
$(1)$ Is the system causal? Why or why not?

I think if $u(t)=0$ for all $t<0$. This means that $u(τ−t)=0$ for all $τ<t$ or, equivalently, for all $t>τ$ and the integrand is zero in range $({-\infty}, t)$.Therefore, we can show that:
$$\int_{-\infty}^\infty x()u(-t) d = \int_{t}^\infty x()u(-t) d =  \int_{t}^\infty x() d$$ So, the system is not causal! Am I right?!

$(2)$ Is the system time-invariant? Why or why not?

Let $ s=t-1$
$$y(s)=\int_{-\infty}^\infty x()u(-s+1) d $$
For invariance, we need to show output resulted from $x(s−s_o)$ = $y(s−s_o)$.
$$y_1(s)=\int_{-\infty}^\infty x(-s_0)u(-s+1) d $$
Now when I change of variable to $z=τ−s_o$, $τ=z+s_0$ , $dz=dτ$ it leads to:
$$y_1(s)=\int_{-\infty}^\infty x(z)u(z+s_0-s+1) dz $$
For $y(s-s_o)$ integral becomes:
$$y_2(s-s_0)=\int_{-\infty}^\infty x()u(-s+1-s_0) d $$
So it seems to time-variant system.
Update:I would appreciate if anyone can answer my question because it's not a homework anymore and yet I didn't find the answer.

Hint
Substitute $s = t-1$. That gets you the equations in a more standard form $y(s) = ...$ Then go through the same excercise.


What does ℜ symbol mean? I came across the symbol $\mathfrak{R}$ in the Kalman introduction: Kalman Filter Intro

$$
x\in \mathfrak{R}^m
$$
It was used in this context:

$$
    x_{k} = Ax_{k–1} + Bu_{k – 1} + w_{k – 1}
$$
What does $\mathfrak{R}$ stand for?

It means $x$ is an $m$-vector of real values. $\mathfrak{R}$ itself is the set of real numbers.


How to achieve Shannon's capacity if coding introduces a spectral efficiency loss? I am a little bit confused concerning different Shannon's notions :

Over an AWGN channel Shannon's capacity is in my knowledge : $R_b<C=B.\log(1+SNR)$.
Is Shannon's "limit" a different notion than that ?

How to convert this equation to a $BER=f(E_b/N_0)$ chart ? for example Here there is a Shannon limit plotted in fig.2, how is it done ?

From the equation stated above one can derive the maximum spectral efficiency : $\eta = R_b/B <\eta_{max} = C/B = \log(1+SNR)$ is that correct ?

If a channel code introduces redundacy bits it introduces a inevitably a loss in $\eta$, so what does Shannon mean when he states that to achieve capacity one has to tend code length to infinity ?  one of my proffessors once said that the LDPC is an example of such a code because it is very long and thus achieves capacity (almost), where is spectral efficiency loss gone ?

You cannot convert the Shannon Capacity formula to a BER "waterfall curve", because the capacity theorem states that in the AWGN channel, there exists a modulation and coding scheme such that the probability of error is arbitrarily small with a capacity bound given by $B \log (1 + SNR)$. Sadly, it doesn't tell us how to construct such a modulation/coding scheme that approaches this bound.

For a code length to approach infinity means that the code depends on an arbitrary number of previous symbols. Heuristically speaking, the "optimal" code will take advantage of correlations in the errors across symbols, so the longer the code is, the more it can exploit this redundancy.


In voice recognition system, how to match digital input to a given word without seconds delay? I've been looking around and found some great C/C++/Java libraries like CMU Sphinx or Windows SAPI. I did not like SAPI, and CMU Sphinx is way too complicated... would need months to understand how it works.

So I would like to start writing my own library, even it would be crap, but to understand all the parts how to recognize a word in real-time.

Now comes in my head a question,... how to match digital input to a given word? How to make it without seconds delay?

So what I understand is that I pass the input into buffer and I have to process it. First I try to create some spectrogram (time, frequency, amplitude) and then to use Fast Fourier Transform. In the database/list I will have stored words with their phonemes.

But how do I match these phonemes with FFT output? I simply do not understand all the steps.

There is a massive range of possible sounds that a human might recognize (or mis-recognize) as the same word in their language.  The problem is searching a large dimensional space representing all those possibilities (continuous variations of spectral and timing differences) for the closest match, perhaps also given some context around the search target for the most probable of multiple near matches.  This is a computational complexity problem that can sometimes be reduced in time by more computational power (a 1000 processor server or supercomputer might be able to do it in less than "seconds" of delay).

Or you can search over less degrees-of-freedom with less robustness or accuracy, and perhaps miss most possible word matches.


Why is the difference of gaussians scale space scale invariant? I'll use the Scale-invariant feature transform algorithm as an example here.  SIFT creates a scale space based on scaled gaussian filtering of an image, and then computes the difference of gaussians to detect potential interest points.  These points are defined as the local minima and maxima across the difference of gaussians.

It is claimed that this approach is scale invariant (among other puzzling invariances).  Why is this?  It's unclear to me why this is the case.

The term "scale-invariant" means the following here.
Let's say you have image I, and you have detected a feature (aka an interest point) f at some location (x,y) and at some scale level s.  Now let's say you have an image I', which is a scaled version of I (downsampled, for instance).  Then, if your feature detector is scale-invariant, you should be able to detect the corresponding feature f' in I' at the corresponding location (x',y') and corresponding scale s', where (x, y, s) and (x', y', s') are related by the appropriate scaling transformation.

In other words, if your scale-invariant detector has detected a feature point corresponding to someone's face, and then you zoom in or out with your camera on the same scene, you should still detect a feature point on that face.

Of course, you would also want a "feature descriptor" which would allow you to match the two features, which is exactly what SIFT gives you.

So, at the risk of confusing you further, there are two things that are scale-invariant here.  One is the DoG interest point detector, which is scale-invariant, because it detects a particular type of image features (blobs) irrespective of their scale.  In other words, the DoG detector detects blobs of any size.  The other scale-invariant thing is the feature descriptor, which is a histogram of gradient orientation, which stays more or less similar for the same image feature despite a change in scale.

By the way, the difference of Gaussians is used here as an approximation to the Laplacian-of-Gaussians filter.


What is the difference between phase delay and group delay? I'm studying some DSP and I'm having trouble understanding the difference between phase delay and group delay.

It seems to me that they both measure the delay time of sinusoids passed through a filter.

Am I correct in thinking this?
If so, how do the two measurements differ?
Could someone give an example of a situation in which one measurement would be more useful than the other?

UPDATE

Reading ahead in Julius Smith's Introduction to Digital Filters, I've found a situation where the two measurements at least give different results: affine-phase filters. That's a partial answer to my question, I guess.

First of all the definitions are different:

Phase delay: (the negative of) Phase divided by frequency
Group delay: (the negative of) First derivative of phase vs frequency

In words that means:

Phase delay: Phase angle at this point in frequency
Group delay: Rate of change of the phase around this point in frequency.

When to use one or the other really depends on your application. The classical application for group delay is modulated sine waves, for example AM radio. The time that it takes for the modulation signal to get through the system is given by the group delay not by the phase delay. Another audio example could be a kick drum: This is mostly a modulated sine wave so if you want to determine how much the kick drum will be delayed (and potentially smeared out in time) the group delay is the way to look at it.


Why do we need to have unity gain at DC for a lowpass filter? I have seen many times, including in a post on this forum, that it is recommended to have unity gain at DC for a lowpass filter. I would like to now the mathematical explanation.

There is no mathematical explanation and it's typically just a convenience for most applications. If you scale so that it's unity at 0 Hz, unity is also the maximum gain at any other frequency so you are less likely to clip or overdrive something in the signal chain (it's still possible though).

It's often based on the assumption or requirement that "frequencies below the cutoff should be unaffected" by a low pass filter. The closest way to get this is to have unity gain at DC.


Response of a system using DFT Let's say that we have two sequences, input sequence $x(n) = [0121]$ and impulse response of a given system $h(n) = [0, 1, -1, 1]$. I need to find response of this system to given input sequence. After that, i need to calculate linear convolution of given sequences.

If we denote the response as $y(n)$ we have $y(n)=h(n)*x(n)$, which means that, due to convolution theorem, in frequency domain we have $Y(k)=H(k)X(k)$. From this, we can find $y(n) = IDFT(Y(k))$.

Considering the fact that i need to find convolution of given sequences, meaning $y(n)=h(n)*x(n)$, that would mean that it should yield same result as when i was doing this by using DFT. However, my final results don't match at all.

$x(n) = [0121] \Rightarrow X(k) = [4, -2, 0, -2] \\ h(n) = [0, 1, -1, 1] \Rightarrow H(k) = [1, 1, -3, 1] \\ Y(k)=X(k)H(k) = [4, -2, 0, -2] \Rightarrow y(n)=[0, 1, 2, 1]$

On the other hand, convolution of given sequences gives the following result:

$y(n) = h(n)*x(n) = [0, 0, 1, 1, 0, 1, 1]$

Not only that result is completely wrong, but also dimensions of vectors i got as a result are not the same. What am i doing wrong? Any help appreciated!

It will match for circular convolution modulo $N$, where $N$ is 4 here. For finite length sequences product of DFT of 2 sequences is equivalent to DFT of circular-convolution of the 2 sequences.

>> cconv([0,1,2,1], [0,1,-1,1], 4)

ans =

 0     1     2     1

>>ifft(fft([0,1,2,1]).*fft([0,1,-1,1]))

ans =

 0     1     2     1


Prefactors computing PSD of a signal with numpy.fft VS. scipy.signal.welch The power spectral density St of a signal u may be computed as the product of the FFT of the signal, u_fft with its complex conjugate u_fft_c. In Python, this would be written as:

import numpy as np

u = # Some numpy array containing signal
u_fft = np.fft.rfft(u-np.nanmean(u))

St = np.multiply(u_fft, np.conj(u_fft))

However, the FFT definition in Numpy requires the multiplication of the result with a factor of 1/N, where N=u.size in order to have an energetically consistent transformation between u and its FFT. This leads to the corrected definition of the PSD using numpy's fft:

St = np.multiply(u_fft, np.conj(u_fft))
St = np.divide(St, u.size)

On the other hand, Scipy's function signal.welch computes the PSD directly from input u:

from spicy.signal import welch

freqs_st, St_welch = welch(u-np.nanmean(u),
     return_onesided=True, nperseg=seg_size, axis=0)

The resulting PSD, St_welch, is obtained by performing several FFTs in segments of the array u with size seg_size. Thus, my question is:

Should St_welch be multiplied by a factor of 1/seg_size to give an energetically consistent PSD? Should it be multiplied by 1/N? Should it not be multiplied at all?

PD: Comparison by performing both operations on a signal is not straightforward, since the Welch method also introduces smoothing of the signal and changes the display in the frequency domain.

Information on the necessity of the prefactor when using numpy.fft :

Journal article on the matter

I would like to post here the link to the accepted answer that I obtained to the question in StackOverflow, as I think it is relevant to this forum.

The original answer by francis is here: post in StackOverflow


MUSIC algorithm derivation Setup

Suppose we have a complex $L\times 1$ signal $\mathbf{x}$ with two tones at (unknown) frequencies and phases defined as:
$$
x_n = A_1 e^{j \omega_1n + \varphi_1} + A_2 e^{j \omega_2n + \varphi_2}
$$
for $0 \leq n \leq L-1$. We observe $\mathbf{x}$ corrupted by additive white Gaussian noise:
$$
\mathbf{y} = \mathbf{x} + \mathbf{w}
$$
where $\mathbf{w} \sim \mathcal{CN}(0, \sigma^2 \mathbf{I})$.

Problem

I'm trying to show a connection between the signal amplitudes and eigenvalues from the eigendecomposition
of the autocovariance matrix $\mathbf{R_y}$ of $\mathbf{y}$ (this is related to MUltiple SIgnal Classification i.e. MUSIC algorithm which splits the $L$ dimensional space into a 2 dimensional signal subspace and an $L-2$ dimensional noise subspace).

$$
\mathbf{R_y} := \mathbf{E}[\mathbf{y}\mathbf{y}^H] = \mathbf{x}\mathbf{x}^H + \sigma^2 \mathbf{I}.
$$

Let $\mathbf{s}_1 = \begin{bmatrix} 1 & e^{j\omega_1} & e^{j2\omega_1} & \ldots & e^{j(L-1)\omega_1}\end{bmatrix}^T$
and similarly define $\mathbf{s}_2$. Then,
$$
\mathbf{x} = \begin{bmatrix}
           \mathbf{s}_1 & \mathbf{s}_2
         \end{bmatrix}
         \begin{bmatrix}
         A_1 e^{j\varphi_1} \\
         A_2 e^{j\varphi_2}
         \end{bmatrix}
$$
and the signal autocovariance matrix can be written as:
\begin{eqnarray*}
\mathbf{x}\mathbf{x}^H &=& \begin{bmatrix}
           \mathbf{s}_1 & \mathbf{s}_2
         \end{bmatrix}
         \begin{bmatrix}
         A_1 e^{j\varphi_1} \\
         A_2 e^{j\varphi_2}
         \end{bmatrix}
         \begin{bmatrix}
         A_1 e^{-j\varphi_1} & A_2 e^{-j\varphi_2}
         \end{bmatrix}
         \begin{bmatrix}
           \mathbf{s}_1^H \\
           \mathbf{s}_2^H
         \end{bmatrix} \\
         &=& \underbrace{\begin{bmatrix}
           \mathbf{s}_1 & \mathbf{s}_2
         \end{bmatrix}}_{\mathbf{S}}
         \underbrace{\begin{bmatrix}
         A_1^2 & A_1 A_2 e^{j(\varphi_1-\varphi_2)} \\
         A_1 A_2 e^{-j(\varphi_1-\varphi_2)} & A_2^2
  \end{bmatrix}}_{\mathbf{M}}
  \underbrace{\begin{bmatrix}
           \mathbf{s}_1^H \\
           \mathbf{s}_2^H
         \end{bmatrix}}_{\mathbf{S}^H}
\end{eqnarray*}
where we note that $\mathbf{M}$ is Hermitian symmetric and therefore has a decomposition $\mathbf{M} = \mathbf{Q\Lambda Q}^H$ where $\mathbf{\Lambda}$ is
a 2x2 diagonal matrix and $\mathbf{Q}^H\mathbf{Q} = \mathbf{Q}\mathbf{Q}^H = \mathbf{I}$.

Therefore,
\begin{eqnarray*}
\mathbf{R_y} &=& \mathbf{SQ \Lambda Q^HS^H} + \sigma^2 \mathbf{I}.
\end{eqnarray*}

Now I am stuck. I suspect I should be able to eventually express $\mathbf{R_y} = \mathbf{PDP^H}$ where the diagonal matrix $\mathbf{D}$ should have two large eigenvalues related to the signal amplitudes $A_1$ and $A_2$, but since $\mathbf{SS^H}\neq \mathbf{I}$ I have no way to absorb the $\sigma^2$'s into the diagonal matrix $\mathbf{\Lambda}$.

The piece I was missing was the distribution of the initial phase values $\varphi_1$ and $\varphi_2$. It is standard to assume that these are uniformly distributed [^]. This leads to:
$$
\mathbf{R_x} = \mathbf{E[xx^H]} = \mathbf{S \Lambda S^H}
$$
where
$\mathbf{\Lambda} = \begin{bmatrix} A_1^2 & 0 \\ 0 & A_2^2
\end{bmatrix}$
and $\mathbf{S} = \begin{bmatrix}\mathbf s_1 & \mathbf s_2 \end{bmatrix}$ so that
$$
\mathbf{R_y} = \mathbf{S \Lambda S^H} + \sigma^2 \mathbf{I}.
$$

Now note that $\mathbf{R_y}$ is Hermitian symmetric and must possess an eigendecomposition with eigenvalues $\lambda_1 \geq \lambda_2 \geq \lambda_3 \geq \cdots \geq \lambda_L$. However, also note that $\mathbf{R_x}$ has two nonzero eigenvalues, say, $\tilde\lambda_1 \geq \tilde\lambda_2$. Since $\mathbf{R_y} = \mathbf{R_x} + \sigma^2 I$ it is easy to verify that $\lambda_1 = \tilde \lambda_1 + \sigma^2$, $\lambda_2 = \tilde \lambda_2 + \sigma^2$ and $\lambda_i = \sigma^2$ for $3\leq i \leq L$.

In the eigendecomposition of $\mathbf{R_y}$, let $\mathbf e_1$ and $\mathbf e_2$ be the eigenvectors corresponding to $\lambda_1$ and $\lambda_2$ and $\mathbf g_1, \ldots, \mathbf g_{L-2}$ be the eigenvectors corresponding to the eigenvalues $\lambda_3, \lambda_4,\ldots, \lambda_L$. One can now show that $\mbox{span}\{\mathbf e_1, \mathbf e_2\} = \mbox{span}\{\mathbf s_1, \mathbf s_2\}$. This is often called the "signal subspace" and the orthogonal subspace $\mbox{span}\{ \mathbf g_1,\ldots \mathbf g_{L-2}\}$ is called the "noise subspace."

Here's the punchline: To estimate $\omega_1$ and $\omega_2$ we need to solve for $\omega$ in $\mathbf s(\omega)^H \mathbf G \mathbf G^H \mathbf s(\omega) = 0$, where $\mathbf G = \begin{bmatrix} \mathbf g_1 & \ldots & \mathbf g_{L-2}\end{bmatrix}$, and $\mathbf s(\omega) = \begin{bmatrix} 1 & e^{j\omega} & e^{j2\omega} & \ldots & e^{j(L-1)\omega} \end{bmatrix}$.

Reference

[^] P. Stoica, R. Moses, Spectral Analysis of Signals 1st Ed. Ch. 4.


Why is my bode plot showing periodic behavior? I found something very interersting when I was trying to simulate a discrete time system. So first I discretized the continuous system $G(s)$ using "c2d" function and "dbode". The resulting bode plot for the discretized system is displaying some sort of periodic behaviour after reaching about 3 hz. So is the phase plot.

Can someone explain why this maybe happening?

When you are discretizing you make transformation $s \to e^{s T}$. This transforms the left half plane into unit disc. But actually it only transforms the $(-j\omega_s/2, j\omega_s/2)$ band and repeats itself by circling for outside of the band, where $\omega_s = 2 \pi / T$. This is because of the periodicity of exponential function.

When you plot the frequency response, you plot $|G(j\omega)|=|\hat{G}(e^{j\omega T})|$ for all $\omega \in [0, \infty)$ where $\hat{G}$ is the discrete transfer function. Therefore the frequency response of any sampled system is periodic.


What is the effect of quantization bucket width in Locality Sensitive Hashing? In the LSH Tutorial by Slaney, it is written that

  Increasing the quantization bucket width w will increase the number of
  points that fall into each bucket. To obtain our final nearest
  neighbor result we will have to perform a linear search through all
  the points that fall into the same bucket as the query, so varying w
  effects a trade-off between a larger table with a smaller final linear
  search, or a more compact table with more points to consider in the
  final search.

As per my understanding, if the bucket size is large, then there should be less number of tables (each bucket is a table right ? ). Since, each point is hashed to a cuket. If a query is hashed to a bucket, then with a linear scan on a larger table, wouldn't it encounter more points than a compact table ?

I think you've got the point, but that you are getting confused by the additional tables that you think should be used to implement buckets.

Let's rephrase the quote from the tutorial first:

  In LSH, you create a partition of the entire space. If you use a large
  bucket width $w$, then obviously you need less buckets to span the
  entire space and thus you get shorter codes. Since you have less
  buckets, your hash table has less slots (obvious) but there are also
  more points per bucket.

Note that by definition of LSH all the points $x$ that fall into a given bucket $B$ have the same binary code:
$$\forall x \in B, \text{LSH}(x) = b,$$
where $b$ is some constant.

Now, let's consider what happens when you make a query on the hash table: if you use $b$ as query to the hash table, you will get all the points $x$ in the above formula.
The only information that you have on these points is that their binary code is equal to $b$, you do not have additional information: they may be all around the bucket, or in its center, etc. So, you can't avoid comparing your exact query point with each of them (using their real coordinates) if you want to find a special neighbour (usually the closest one).

Since you need to loop over all the points inside the bucket, this search is linear, and here comes the trade-off mentioned in the tutorial.

When you implement some LSH, you are free to use whatever technique you want for both the hash table and the buckets.
In my hashing implementation, the hash table was actually a table (for simplicity, because finding a bucket becomes then a simple table lookup), but the buckets themselves were implemented as lists for memory efficiency (you allocate a slot only when you add a point to the bucket).

Your question implies that you would implement the buckets themselves as arrays, but this may not always be desirable when you do not know anything about the distribution of the points in your data set or when your programming framework (in my case, vanilla C/C++) does not allow you to easily resize arrays.


fft result order compared to direct implementation of DTFS using text book To find the complex Fourier coefficients of this sequence, then The $k^{th}$
coeff. $c[k]$ is given by the DTFS formula

$$
c[k]=%
%TCIMACRO{\dsum \limits_{n=0}^{N-1}}%
%BeginExpansion
{\displaystyle\sum\limits_{n=0}^{N-1}}
%EndExpansion
x[n]\exp(-in\frac{2\pi}{N}k)
$$

(Lets forget about the $1/N$ scaling for now). My question is on the ordering
of the result.

For example, given $x=[1,1,0,0]$, hence $N=4$. Then applying the above gives

$c[0]=2,c[1]=1-i,c[2]=0,c[3]=1+i$

but Matlab's  fft  gives

c=fft([1 1 0 0]);

$c[0]=2,c[1]=1+i,c[2]=0,c[3]=1-i$

So matlab's $c[1]$ is the formulas $c[3]$ and Matlab's $c[3]$ is the formulas
$c[1]$. So the order seems to be swapped around.

Ofcourse they are still complex conjugate of each others. So both answers are correct.

My question is on the order of terms and why they come out this way.

Any thought on what that is, and why the result is not in the same order as
the formula gives? And if it is possible to make fft give the same order as the formula?

add

I found the problem. I used transpose to change row result from fft to column to make it easy to see. But I did not use the correct tranpose.

When I transposed the result, which was a row, to a column to easier to see, the terms got switch around. Here is an example

z=[1 2+1i 3 4+1i];
z'

   1.0000 + 0.0000i
   2.0000 - 1.0000i
   3.0000 + 0.0000i
   4.0000 - 1.0000i

I should have used non-conjugate tranpose:

z=[1 2+1i 3 4+1i];
z.'

   1.0000 + 0.0000i
   2.0000 + 1.0000i
   3.0000 + 0.0000i
   4.0000 + 1.0000i

So, terms are in correct order all along.

User error. sorry.

Octave returns your same first result:

octave-3.2.3:1> fft([1 1 0 0])
ans =

   2 + 0i   1 - 1i   0 + 0i   1 + 1i

While 4*ifft([1 1 0 0]) results in what you claim is Matlab's result.

Please check you are not confusing direct and inverse FFT.


Normalized Frequency to rads/sec, Hz and discrete/digital frequency I have a signal sampled at 16000 Hz
I need to convert the normalized frequency of Matlab 0.325 to radians/sec Hz and Discrete frequency.
Can anyone explain this to me?
And btw Matlab plots the normalized frequency as $\pi$ rads/sample. So we have to divide the normalized frequency with $\pi$?

The continuous-time frequency $\Omega$ and "discrete" frequency $\omega$ given a sampling period of $T_s$ is given by
$$\omega = \Omega T_s = 2{\pi}fT_s \space\text{rads/sample}$$
Where $-\pi \lt\omega\le\pi $ assuming double-sided spectrums.
It's important to note the units: the continuous frequency $\Omega$ is in rads/s, so multiplying by the sampling period gives you the unit of rads/sample.
If you rewrite the above expression using the sampling frequency $f_s$ you get
$$\omega = 2{\pi}\frac{f}{f_s} \space\text{rads/sample}$$
MATLAB divides this by $\pi$ so that we get a new set of normalized frequencies $\omega'$
$$\omega' = \frac{\omega}{\pi} =  2\frac{f}{f_s}$$
When you sample at a rate of $f_s$, your single-sided unaliased spectrum would be in the range of $[0, \frac{f_s}{2}]$. If you take this range and use the equations for $\omega$ and $\omega'$ you get the ranges

$[0, \pi]$ for $\omega$
$[0, 1]$ for $\omega'$

The second range is what MATLAB uses as the default when plotting discrete frequency responses. In order to report the frequency as "rads/sample" the factor of $\pi$ is re-introduced which is why you see the axis labeled as "x $\pi$ rads/sample". This convention is used to help you quickly identify frequency values as related to your sampling frequency $f_s$ as we'll see in some examples blow.
As an example, a normalized discrete frequency of 0.5 from MATLAB at your sample rate of 16 kHz gives you the continuous time frequency of
$$f = \frac{\omega'f_s}{2} = \frac{(0.5)16000}{2} = 4 \text{ kHz}$$
A value of 1 will yield
$$f = \frac{\omega'f_s}{2} = \frac{(1)16000}{2} = 8 \text{ kHz}$$
The latter is exactly the Nyquist frequency of your signal and gives us the analog to discrete domain frequency mappings we expect. You can find more information here.


Sample time in convolution I have a signal sampled at a period of 0.01 and a channel impulse response sampled at 0.1250. I convolved the signal with the  channel impulse response.

It is right to perform the convolution between two signals with diffrent sample times?
If it is,  how can I calculate the sample time of the convolution result?

It is right to perform the convolution between two signals with diffrent sample times?

No, aside from a very few special use cases (if you have to ask, not yours), you need to resample either of them so that both are at the same sample rate.

Also notice that sampling a channel impulse response at 8 Hz only gives you information about these 8 Hz (or 4 Hz, depending on whether this is real direct sampling or a sampled equivalent complex baseband representation of the channel). You're trying to apply it to a 100 Hz (or 200 Hz, depending...) wide channel: that won't give you anything useful, even if you resampled. Your channel impulse response measurement simply isn't sufficient.


How to calculate FFT of a non-uniformly sampled signal? Is it possible to calculate the FFT of a signal that was not sampled uniformly? I have an MPU6050 connected to an Arduino uno.
I write the samples at the default frequency, approximately 20 milliseconds,
but it's not always constant (see the image). Also see this related question.

This is my first time I approach to the FFT, and I have a civil engineer background, so is my first time with signal analysis, never studied Signal Theory before.

My interest for the FFT is to define the best low pass filter (example apply a butterworth but I don’t know how to choose the filter order and cutoff frequency.

Thanks for your time, thanks for you patience. I tried to solve this with the code below, but I think is only rubbish.

M=csvread('DATALOG.TXT');
n=length(M);
cleaning=500;
n=(n-cleaning);

MatFilt= M(cleaning:n,:);
m=length(MatFilt);
A= MatFilt(2:m,1);
B=MatFilt(1:m-1,1);
C=A-B;
meanrate= mean(C);
f= 1000/meanrate;
T= sum(MatFilt);
TemporalSample= MatFilt(:,1);
Ax=MatFilt(:,2);
Axs=Ax*19.6/32768; %I converted the data in m/s^2
SignalLength= (MatFilt(m,1)-MatFilt(1,1))/1000;
plot(TemporalSample,Axs);
pause
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('time (milliseconds)');
NFFT = 2^nextpow2(m); % Next power of 2 from length of y
Y = fft(Axs,NFFT)/m;
f = f/2*linspace(0,1,NFFT/2+1);

You are essentially referring to nonuniform signal sampling in which data samples are not acquired at exact periods but have same random jitter deviations from their exact uniform timing positions.

Now under some suitable conditions (such as the jitter length being less than a period), there are very effective algorithms which can perfectly convert the nonuniformly sampled data into their equivalent uniform samples. (i.e., given the values of nonuniform samples they give a new set of values which correspond to the perfect uniform samples).

After this nonuiform to uniform conversion, you can then apply the usual FFT, which assumes that the data to be transformed was uniformly sampled.

However, before such a DSP attempt, I would suggest you understand the sampling and read out mechanism of the Arduino Uno and MPU 6050 6DOF sensor combination. That sensor can be read in either a single data read mode or in a burst buffer data read mode. In this second mode, I think it will provide a set of buffered data which are sampled by the MPU6050 unit itself as uniform as possible, compared to a sigle data read mode mastered by the Arduino.

From your data, there seems to be a drift of a few milliseconds per sample. I assume a 20 milli second of sampling period is quite achievable with less than 1 ms of timing accuracy using the Arduino in a proper configuration. So may be your main problem is on the software configuration?

EDIT: Try the following MATLAB code as is:
NOTE: The following code MAY result in WORSE results than available before processing... The success highly depends on the application.

% NONUNIFORM SIGNAL DEFINITION
% ----------------------------
N=100;       % Converts a block of length N nonuniform samples to uniform
K=N;        %
M = N*K;    %
Tsam = 0.020; % Sampling period in seconds,

% tp represents the nonuniform sampling times
tp = [0:Tsam:(K-1)*Tsam] + Tsam*(rand(1,K)-0.5);
xnunds = sin(2*pi*1*tp);        % A simulated nonuniformly sampled data...

% GENERATE THE MATRIX AND ITS INVERSE
% -----------------------------------
A=zeros(K,K);       % initial matrix to be inverted
for p=1:K
    for q=1:K

        A(p,q)=sinc((tp(p)-tp(q))/Tsam);

    end
end
AI=inv(A);          % Inverse matrix is here

% APPLY THE CONVERSION FORMULA
% -----------------------------
xint=zeros(1,K);    % interpolate to uniform samples

for n=1:N:M

    ACS=0;

    for m=1:K       % per nonuniform sample loop

        PSI=0;

        for q=1:K   % calculation of PSI for each m

            PSI = PSI + AI(q,m) * sinc( ((n-1)*Tsam/N - tp(q))/Tsam);

        end%q

        ACS=ACS+xnunds(m)*PSI;  % acumulation of multiplicant

    end%m

    xint( (n-1)/N + 1)=ACS;    % result for each n is assigned
end%n

figure,stem(xnunds);
hold on
stem(sin(2*pi*1*[0:K-1]*Tsam),'g');
title('nonuniform samples imposed on true samples');

figure,stem(xint,'r');
hold on;
stem(sin(2*pi*1*[0:K-1]*Tsam),'g');
title('interpolated uniform samples imposed on true uniform samples)


Creating a digital filter, from Laplace to $\mathcal Z$-transform (zero order hold) to code? I'm trying to create a digital filter in code(C) but any language is fine. Now I've got an analogue filter that I have represented by an equation in the Laplace domain and I want to try and implement it digitally.

So my filter has this form in the Laplace domain:
$$\frac{as+b}{cs^2+ds}$$

I then use MATLAB's c2d command which uses the zero order hold transformation (I have a really poor grasp on this, so this might be wrong) and it gives me this formula:

$$\frac{\left(5\cdot 10^5\right)z-67}{z^2-z}$$

I tried following an example that I found that used the Tustin's method, though when I use the c2d function in MATLAB with Tustin it gives me an error.

My attempt has been

$$\frac{hz-i}{jz^2-kz}$$

$b_0=-i, b_1=h, b_2=0, a_0=0, a_1=-k, a_2=j$

Then from this I've tried (which is wrong)
\begin{align}
\text{output}&=z_0 b_0+z_1b_1+z_2b_2\\
z_2&=z_1\\
z_1&=z_0\\
z_0&=\text{input}-a_0z_0-a_1z_1-a_2z_2
\end{align}

The example I looked at used a tustin or bilinear conversion not a zero order hold(the default for matlabs "c2d" command). So this is more an answer to what i wanted to do rather than the question that i asked above.

I solved the following (converting the s domain function into code) by taking the s domain function.
$$\frac{as+b}{cs^2+ds}$$

and putting this into matlab (matlab command "g=tf([a b],[c d 0])"). Then performed the bilenear conversion with the matlab command "c2d(g,Ts,'tustin')" where g was my transfer funtion and Ts my sampling rate. This produced the output

$$\frac{ez^2+fz+g}{iz^2+jz+k}$$

The a and b coeficients can then be taken from this equation such that(if $i!=1$ the equation needs to be multiplied through by the inverse of "i"):
$b0=e$ $b1=f$  $b2=g$
 $a0=i$  $a1=j$ $a2=k$

this can then be converted to code by setting the initial states for simplicity let $$z0=z1=z2=0$$

then set up a loop that repeats the following algorithm

$$output=z0*b0+z1*b1+z2*b2$$
$$z2=z1$$
$$z1=z0$$
$$z0=input-a1*z1-a2*z2$$

For anyone else that got lost like me, this is known as an IIR filter and googling IIR filter design helped sooo much.


Do you recognize this discrete realization of a low pass filter? I have a simple 1st order lag filter implementation that performs well, shown in  the block diagram below.  This is legacy 'code' we've been using for many years, the original developers have long since moved on.  I've been asked to doc the requirements (development of the difference equation)  starting from the s domain and assuming $\frac{Y(s)}{U(s)} = \frac{1}{\tau s+1}$.
$k1$ is $0.5+\frac{\tau}{T_s}$, $k2$ is $0.5$. $\tau$ is the time constant, in sec, and $T_s$ is the time between samples, in sec.
My assumption was I would arrive at a bilinear transform of $\frac{1}{s+1}$ but it's just not happening as the topology is wrong.  Any hints?
Thanks,
Chris

This is a discrete time filter and you have a block diagram, so forget the s-domain, and just work directly with difference equations and the Z transform.
In your block diagram, let's call the signal going into the delay block $w[n]$ and the signal coming out of the division block $v[n]$.
Writing some initial difference equations:
$$\begin{align*}
w[n] &= w[n-1] + v[n] \quad \Rightarrow \quad v[n] = w[n] - w[n-1]\\
\\
v[n] &= \dfrac{1}{k_1} u[n] - \dfrac{1}{k_1} w[n-1] \\
\\
y[n] &= w[n] - k_2 v[n]\\
\end{align*}$$
Performing some simple algebra, we get:
$$\begin{align*}
y[n] &= \left(1 - k_2\right) w[n] +k_2 w[n-1]\\
\\
w[n] - \left(1 - \dfrac{1}{k_1}\right)w[n-1] &= \dfrac{1}{k_1}u[n]\\
\end{align*}$$
So take the Z transforms and manipulate:
$$\begin{align*}
Y(z) &= \left(1 - k_2\right) W(z) +k_2 z^{-1}W(z)\\
\\
\dfrac{Y(z)}{W(z)} &= \left(1 - k_2\right)  +k_2 z^{-1} \\
\\
\dfrac{Y(z)}{W(z)} &= \left(1 - k_2\right) \left(1 - \dfrac{k_2}{k_2-1}z^{-1}\right)\\
\\
W(z) - \left(1-\dfrac{1}{k_1}\right)z^{-1} W(z) &= \dfrac{1}{k_1}U(z)\\
\\
\dfrac{W(z)}{U(z)} &= \dfrac{1}{k_1} \cdot \dfrac{1}{1-\dfrac{k_1 -1}{k_1}z^{-1} }\\
\end{align*}$$
The final transfer function is then:
$$\begin{align*}
\dfrac{Y(z)}{U(z)} &= \dfrac{Y(z)}{W(z)}\cdot\dfrac{W(z)}{U(z)}\\
\\
&= \dfrac{1-k_2}{k_1} \cdot \dfrac{1 - \dfrac{k_2}{k_2-1}z^{-1}}{1-\dfrac{k_1 -1}{k_1}z^{-1} }\\
\\
&= \dfrac{0.5}{0.5 + \frac{\tau}{T_s}}\cdot \dfrac{1 - (-1)z^{-1}}{1-\dfrac{\left(0.5 + \frac{\tau}{T_s}\right) -1}{0.5 + \frac{\tau}{T_s}}z^{-1} }\\
\end{align*}$$
You have a zero at the Nyquist frequency, $-1$ on the unit circle, so this is definitely looking like a lowpass filter.
You have a pole on the real axis just inside the unit circle, near $1$, at
$$ z_p = \dfrac{\left(0.5 + \frac{\tau}{T_s}\right) -1}{0.5 + \frac{\tau}{T_s}} = 1 - \dfrac{1}{0.5 + \frac{\tau}{T_s}} $$
which should provide a strong response near DC.  This is what you would expect from a lowpass filter.
You can do some algebra and take an inverse Z-transform to get the final difference equation:
$$\begin{align*}
\left( 1-\dfrac{\left(0.5 + \frac{\tau}{T_s}\right) -1}{0.5 + \frac{\tau}{T_s}}z^{-1}\right)Y(z) = \dfrac{0.5}{0.5 + \frac{\tau}{T_s}}\cdot \left(1 - (-1)z^{-1} \right) U(z) \\
\\
y[n] - \dfrac{\left(0.5 + \frac{\tau}{T_s}\right) -1}{0.5 + \frac{\tau}{T_s}} y[n-1] = \dfrac{0.5}{0.5 + \frac{\tau}{T_s}} u[n] + \dfrac{0.5}{0.5 + \frac{\tau}{T_s}} u[n-1]\\
\\
y[n] = \left(1-\dfrac{1}{0.5 + \frac{\tau}{T_s}}\right) y[n-1] + \dfrac{0.5}{0.5 + \frac{\tau}{T_s}} u[n] + \dfrac{0.5}{0.5 + \frac{\tau}{T_s}} u[n-1] \\
\end{align*}$$
So the filter is a scaled average of the two most recent input values with decaying feedback of the output value.
This could have been implemented with a much more straightforward block topology than the wonky one you were given.
In fact this difference equation is a 2 sample Moving Average filter, followed by an Exponentially Weighted Moving Average (EWMA) filter:
$$\begin{align*}
x[n] &= 0.5u[n] + 0.5u[n-1] \quad \text{(Moving Average)}\\
\\
y[n] &= \left(1-\dfrac{1}{0.5 + \frac{\tau}{T_s}}\right) y[n-1] + \dfrac{1}{0.5 + \frac{\tau}{T_s}} x[n] \quad \text{(EWMA)} \\
\end{align*}$$
From this answer: https://dsp.stackexchange.com/a/40465/28112 the -3 dB cutoff frequency of your EWMA filter is:
$$f_{3dB} = \dfrac{1}{2\pi T_s} \cos^{-1}\left[\dfrac{\left(\dfrac{1}{0.5 + \frac{\tau}{T_s}}\right)^2}{2\left(1-\dfrac{1}{0.5 + \frac{\tau}{T_s}}\right)}\right]$$
Two other observations for you:
$k_2$ sets the balance of the weight given to the inputs $u[n]$ vs. $u[n-1]$
$k_1$ sets the $\alpha$ of the EWMA filter, which sets the balance of the weight given to the previous output vs. the new input.
Update
So now that we know where the digital filter's pole is, let's work through from the continuous time filter's pole, to see if it matches, using the bilinear transform relationship.
$$\begin{align*}
s_p = -\dfrac{1}{\tau} &= \dfrac{2}{T_s} \cdot \dfrac{z_p -1}{z_p+1}\\
\\
z_p+1 &= -2\dfrac{\tau}{T_s}\left(z_p-1\right) \\
\\
z_p\left(1 + 2\dfrac{\tau}{T_s}\right) &= -1 + 2\dfrac{\tau}{T_s} \\
\\
z_p\left(1 + 2\dfrac{\tau}{T_s}\right) &= -2 + \left(1 + 2\dfrac{\tau}{T_s}\right) \\
\\
z_p\left(0.5 + \dfrac{\tau}{T_s}\right) &= \left(0.5 + \dfrac{\tau}{T_s}\right) -1 \\
\\
z_p & = \dfrac{\left(0.5 + \dfrac{\tau}{T_s}\right) -1}{0.5 + \dfrac{\tau}{T_s}} = \dfrac{k_1 -1}{k_1}
\end{align*}$$
So that checks out. Someone did apparently start from a single pole filter of the form
$$H_a(s) = g\cdot\dfrac{1}{1+\tau s}$$
$H_a(s)$ has a zero as $s \rightarrow \infty$ (along any direction from the origin of the s-plane).  Looking again at the bilinear transform relationship:
$$\begin{align*}
s_0 &= \dfrac{2}{T_s} \cdot \dfrac{z_0 -1}{z_0+1}\\
\\
\lim_{z_0 \to -1} s_0 &\rightarrow -\infty\\
\\
z_0 &= -1\\
\end{align*}$$
Which also agrees with the zero we found analyzing the digital filter.
Update 2
And just to give some insight as to why $k_2 = 0.5$
$$\begin{align*}
s_0 &= \dfrac{2}{T_s} \cdot \dfrac{z_0 -1}{z_0 +1}\\
\\
&= \dfrac{2}{T_s} \cdot \dfrac{\frac{k_2}{k_2-1} -1}{\frac{k_2}{k_2-1} +1} \\
&= \dfrac{2}{T_s} \cdot \dfrac{k_2 - k_2 + 1}{k_2 + k_2-1}\\
\\
&= \dfrac{2}{T_s} \cdot \dfrac{ 1}{2k_2 - 1}\\
\\
&= \dfrac{1}{T_s} \cdot \dfrac{ 1}{k_2 - 0.5}\\
\end{align*}$$


Relationship between filters, Fourier magnitude, and frequency I'm just beginning to study Fourier transforms and I'm having trouble getting this concept down.
I know what a low/high pass filter does to a signal.

What is the relationship between Fourier magnitudes (dB) and frequency (Hz)?
How would Fourier magnitudes get affected when a low/high pass filter is passed on a signal?
How about an all pass filter?

(This explanation is going to be rough around the edges so that I can build your intuition. There are subtleties involved, but this should be good enough to start you off).

Forget about the dB scale for a moment. Let us take a step back:

Imagine that you have some signal that occupies a bandwidth from say, 50 Hz to 100 Hz. Its spectrum is flat across this band. That is, it has the same power at 50 Hz, as it does at 51 Hz, as it does at 52 Hz, etc etc. In this example, we say that it has a constant (flat) power-spectral-density within this band. Let us say that the value of the PSD in this band is 1 Watts / Hz. This means that across 1 Hz of frequency, we will measure 1 Watt in power. Across 50 to 100 Hz, we would measure of course, 50 Watts in power. (To compute your PSD values, they are simply the absolute value of your fourier amplitudes, squared).

So far so good.

Now, imagine a filter, that removes frequencies above 75 Hz and also removes frequencies below 50 Hz. It however keeps frequencies between 50 and 75 Hz. To "keep" a frequency, means it multiplies the power of a signal at that frequency by 1. To "remove" a frequency, means it multiplies the power of a signal at that frequency by 0.

So, you have a signal that is 1 Watt/Hz from 50 to 100 Hz. And you pass it through a filter that multiplies all frequency powers from 50 to 75 by 1, ("keeps them"), and it also multiplies all frequency powers less than 50 and above 75 Hz by 0. ("removes them").

If you pass your signal through this filter, what might you get? Well, you just removed all the energy of your signal from 75 to 100 Hz. Previously those frequencies had a power spectral density of 1 Wattage/Hz, however now they are 0. You have just filtered your signal.

If I asked you at this point, "Plot for me, the PSD of your signal after your filtering", you would draw 1's for frequencies in the 50 to 75 Hz range, and draw 0's for all other frequencies. Now if I said, "I dont like this linear scale. I want you to draw this for me on a dB scale. A log scale." This means that you now simply take the log of the PSD axis. So everywhere we had 0 watts/Hz, they are now -infinity watts/Hz dBW/Hz, (db Watts), and everywhere you had 1 Watt/Hz, we now have 0 dBW /Hz.

Taking it so dB scale doesn't change your frequency axis, only your fourier amplitude/power axis.


How can define an indicator which measures the degree of similarity between two signals? The similarity of two signals is calculated by cross correlation. But, how to define an indicator which quantitatively measures the degree of similarity between two signals? Thanks.

The general topic of finding similarities between signals is wide ranging:

are the signals of same sampling, length, offset, shift or scale?
where do they take their values (discrete, real, complex)?
are they stationary? noisy?
what do you consider similar (whole signals, chunks, specific features)?
which are the invariances looked for?
and most important: what is your goal?

The following works provide an overview of some common similarity metrics, which one can normalize into an index:

Clustering of time series data --- a survey, section 2.2. Similarity/distance measures
Similarity Measures and Dimensionality Reduction Techniques for Time Series Data Mining, section 3. Similarity measures


Is $\cos(x[n])$ stable? I am highly confused about the stability of $\cos(x[n])$.

If we provide a bounded input such as $x[n]=u[n]$, the output is bounded.

Now if we provide an unbounded signal $x[n]=\delta[n]$, the output oscillates between $-1$ and $1$. So I think it's marginally stable in this case.

Am I correct? I have this confusion because in my book it says it's unstable when $x[n]=\delta[n]$. So my questions are: Is the function BIBO stable and also is it stable for an unbounded input ?

Note: $u[n]$ and $\delta[n]$ are unit step and dirac delta functions respectively.

Correction: As mentioned in the comment, $x[n]=\delta[n]$ is bounded too.

But this points more strongly to $\cos(x[n])$ being stable, isn't it?

The I/O relation
$$y[n] = \cos(x[n])$$

represents a stable system.

Because by definition of $\cos(x)$, irrespective of the input $x[n]$, its output will always be bounded by $$-1 \leq \cos(x[n]) \leq 1 $$ for all $x[n]$ for all $n$. Even if $x[n]$ is undefined ( such as would be the case with $x[n] = 1/n$ at $n=0$ as Marcus pointed out) we for sure know that the system output is bounded (the value of input is undefined, but the output is known to be within $\pm1$.)


How to implement discretized CWT with FFT? I tried to make circular convolution (correlation actually because it is easier to understand for me) look like the Continuous Wavelet Transform (CWT) but got stuck because the conjugation was on the wavelet and the sign was different too (resembled $x[m](y[m-n])^*$ instead of $x[m] y[n-m]$).

I have a signal sampled from my microphone, and I have numpy and scipy. But since I can't get the indices to play nice I don't know how to implement cwt using fft.

Realized that the (unitary-/ortho-)DFT is basically a list of inner products
$$F(x)_k=\langle x(\_), |n|^{-1/2}\exp(\frac{i\cdot 2\cdot \pi \cdot k \cdot \_}{n})\rangle_k=\langle x(\_), f_k(\_)\rangle_k\\
F^{-1}(X)_k=\langle X(\_), f^*_k(\_)\rangle_k$$
since the second argument gets conjugated in the complex inner product of $\ell^2$.

In a similar fashion I guess a discretized CWT (DCWT) would be something like
$$\Phi(x)_{(a,b)}=\langle x(\_), |a|^{-1/2} \phi(\frac{\_-b}{a}))\rangle_{(a,b)}=\langle x(\_), \phi_{(a,b)}(\_)\rangle_{(a,b)}$$
(haven't found a formula for the inverse that is symmetric like the dft yet so the normalization constant might be wrong).

Convolution would be something like
$$(x\ast y)_{k} = \langle x(\_) , y^*(-\_+k)\rangle_{k}$$
(conjugation of conjugation should be identity. Should it be divided by $n$?)

Cross-correlation would be something like (should it be divided by $n$?)
$$(x\star y)_k = \langle x^*(\_), y^*(\_+k)\rangle_k$$
wolfram also says (I think. Should it be divided by $n$?)
$$
(x\star y)_k = (x^*(-\_)\ast y(\_))_k = \langle x^*(-\_) , y^*(-\_+k)\rangle_{k}
$$

I think the convolution theorem says
$$
[F\cdot F](x,y)_k=F(x)_{k}\cdot F(y)_{k}\\
(x\ast y)_k = F^{-1}([F\cdot F](x,y)_{\_})_k
$$
(should it be multiplied by $|n|^{1/2}$ or $|n|^{3/2}$ when using unitary dft? should it be as is or multiplied by $n$ when using the default dft?)

I think I want the b (or maybe some fraction of a and b) of the DCWT to be the k of the convolution/correlation somehow? Perhaps
$$
y_a(t)=|a|^{-1/2}\phi^*(-a^{-1}\cdot t)\\
\Phi(x)_{(a,b)}=(x\ast y_a)_b
$$

To be honest, if I want to use a CWT, I'll use the toolbox Matlab provides.

I have every confidence that I could code a CWT, but an efficient, debugged, and tested set of routines is not an interest.

There are a number of open source wavelet codes that a Google search will return, such as

https://stackoverflow.com/questions/9606458/looking-for-a-good-c-c-wavelet-library-for-signal-processing

so there is code available that you can use and/or modify or just get ideas from.

Also, you haven't posted what you have so far produced, code or equations.  No one can offer any concrete suggestions or advice given your description.


Preprocessing signal before LPC I have a signal that contains the speaking of three vowels. The signal is contaminated with dc offset and an interference at 60 Hz. I want to perform LPC for the synthesis of the spoken vowels. Do I have to remove the dc offset and the interference before performing LPC?

My initial signal looks like this:

There is also another question of mine about LPC here https://stackoverflow.com/questions/28909899/weird-output-of-lpc-algorithm

Usually before processing the speech signals, they are passed through a high pass filter with cutoff frequency of around 150 Hz to remove such artefacts.


Rectangular 8 QAM bit error rate calculation I am simulating an 8-QAM communication system over an AWGN channel for various signal to noise ratios, in order to plot the average bit error. However I have not been able to derive an equation for the QAM system in terms of the Q-function.

The 8-QAM system is rectangular with 2 possible values for the x basis vector and 4 possible values for the y basis vector such that the constellation consists of 8 symbols:

(-d/2,3d/2), (d/2,3d/2)
(-d/2,d/2), (d/2,d/2)
(-d/2,-d/2), (d/2,-d/2)
(-d/2,-3d/2), (d/2,-3d/2)

where d is the distance between symbols.

The symbols are coded using grey coding for the x and y component as follows:

(0,10), (1,10)
(0,11), (1,11)
(0,01), (1,01)
(0,00), (1,00)

My attempt to derive an expression for bit error rate was based on the symmetry of the constellation as it consists of only two unique error scenarios, since the received vectors are:

r1 = s1 + n1, r2 = s2 + n2

The probability of error in each scenario is
$$
P(bit error scenario 1) = P(r1<0|s1) + P(r2>d|s1) + P(r2<0|s1)
= 3Q(d/2\sigma)
$$
$$
P(bit error scenario 2) = P(r1<0|s1) + P(r2<d|s1)
= 2Q(d/2\sigma)
$$
Average probability of bit error is then:

$P(bit error) = (5/2)Q(d/2\sigma)$

Assuming the transmit signals are equi-probable. and:

$d = \sqrt(2E_s/3), $
$\sigma =\sqrt(N_0/2)$

however there is significant disparity between this theoretical error rate and the simulated error rate. I would appreciate it if someone could point me in the right direction, thanks.

Neglecting square terms of the Q function, your expression is correct, but it is the symbol error probability, not the bit error probability:

$$P_s=2.5\,Q\left(\frac{d}{2\sigma}\right)\tag{1}$$

With Gray encoding, the bit error probability is closely approximated by

$$P_b\approx\frac13Ps\tag{2}$$

because there are $3$ bits per symbol, and the most likely symbol errors (i.e. between neighboring symbols) cause only one bit error.

Furthermore, if $N_0$ is the one-sided noise density, then $\sigma=\sqrt{N_0/2}$, and the energy per bit is $E_b=E_s/3$, where $E_s$ is the symbol energy. It remains to express the distance $d$ between the symbols in terms of $E_s$. The energy of the four inner symbols is $d^2/2$. The energy of the corner symbols is $5d^2/2$. So for equally probable symbols, the average symbol energy is

$$E_s=\frac12\cdot\frac{d^2}{2}+\frac12\cdot\frac{5d^2}{2}=\frac{3d^2}{2}\tag{3}$$

which gives

$$d=\sqrt{\frac{2E_s}{3}}=\sqrt{2E_b}\tag{4}$$

Plugging (4) into (1) gives for the symbol error probablity

$$P_s=2.5\,Q\left(\sqrt{\frac{E_b}{N_0}}\right)\tag{5}$$

which, combined with (2), gives for the average bit error probability

$$E_b\approx\frac{5}{6}\,Q\left(\sqrt{\frac{E_b}{N_0}}\right)\tag{6}$$

Note that Eq. (6) for the bit error probability was derived under the assumption that the SNR is relatively high. This assumption is implicit in Eq. (1), where quadratic terms of the Q function were neglected. The exact expression for the symbol error probability is

$$P_s=2.5\,Q\left(\frac{d}{2\sigma}\right)-1.5\,Q^2\left(\frac{d}{2\sigma}\right)\tag{a}$$

For $d/2\sigma\rightarrow 0$ the expression (a) approaches $7/8$. This is a consequence of all symbols being detected with equal probability, and since only one of them is the correct symbol, the symbol error probability must be $7/8$.

The expression (2) for the bit error probability assumes that each symbol error causes only one bit error, which is true if Gray encoding is used and if it is assumed that symbol errors only occur between neighboring symbols, the latter being a reasonable assumption at relatively high SNRs. However, for $d/2\sigma\rightarrow 0$ any symbol will be confused with any other with equal probability. There is only one symbol (the correct one) where there are zero bit errors. There are $3$ symbols causing one bit error, $3$ symbols causing $2$ bit errors, and one symbol causing $3$ bit errors. Since for very bad SNRs all symbols are detected with equal probability, the average number of wrong bits per detected symbol is

$$\frac18\cdot 0+\frac{3}{8}\cdot 1+\frac{3}{8}\cdot 2+\frac18\cdot 3=\frac{3}{2}$$

And since there are $3$ bits per symbol, the average bit error rate for $d/2\sigma\rightarrow 0$ is

$$E_b=\frac{3}{2}\cdot\frac13=\frac12$$

as expected.


SNR of a pulsed signal in noise I have a digitized acoustic data of length L in which a pulsed signal (chirp) of finite duration and bandwidth is merged starting at some time instant t1 and ending at time instant t2. I wan to calculate SNR of this pulse signal. What is the best way to calculate this SNR?

Thanks & Regards

You detect the chirp, and subtract it from the mixture. The remainder is your signal of interest (plus other sources of noise which you don't mention). Then you calculate the power of each component, and calculate the ratio. That's SNR.

The challenge is usually detecting the chirp, but that can be done, depending on the level of information you have about it (chirp rate, phase of the chirp, beginning of the chirp, bandwidth, amplitude) using simple techniques like a correlation with a prototype chirp, or "blinder" methods like letting a PLL run on the signal and watching when it locks.


Time scaling and shifting of delta function Learning signals and systems.
Solving time scaling and shifting problems.

For the  question
$$x(t) = u(2t - 1)$$
First we shift by 1 to the right side and then we do time scaling , i.e divide by 2 on the time axis.

$x(t) = δ(2t - 1) $

Can we do the same thing for the above impulse function

Is this the correct order of solving this:

Shift to right by 1
time scale by half
change the area of the delta function by multiplying it to 1/2.

Note that the time-scaled and shifted Dirac delta impulse $\delta(2t-1)$ is non-zero where $2t-1=0$ is satisfied, i.e., at $t=\frac12$. It's area is indeed $\frac12$ because we have

$$\delta(at)=\frac{1}{|a|}\delta(t)\tag{1}$$

So for your example you get

$$\delta(2t-1)=\frac12\delta\left(t-\frac12\right)\tag{2}$$