How to fix frequency drift when producing an arbitrary length signal by repeatedly looping over a single-period waveform? I have an audio application that produces an arbitrary length signal by repeatedly looping over a single-period waveform. In pseudo-code:

// Build a single period of a waveform at the given frequency
samples_per_waveform = int(sampling_rate / frequency)
waveform = make_sin_waveform(samples_per_waveform)

...

// Repeat the single period waveform to produce a continuous signal
count = int(duration * sampling_frequency / samples_per_waveform)
repeat `count` times:
     write(waveform)

Since the number of samples per waveform is obviously an integer, I think this is subject to rounding error:

As an example, when rounding to the nearest integer, at a 5kHz sampling rate, a single-period waveform for a 300Hz signal would require 17 samples--and thus leads to an actual frequency of 294Hz. I understand the higher will be the sampling rate, the lower will be the rounding error. However, it will still exist.

Is there a way to fix that, or is this approach doomed from the start?

EDIT: As a possible workaround, I envision to keep track of the rounding error in the output loop, and duplicating one sample when the error exceeds a sample duration. I don't know if this is clear enough, but it is somewhat similar in my mind to what the Bresenham line drawing algorithm is doing.

My suggestion is to use a Numerically Controlled Oscillator (NCO). A typical NCO effectively has one cycle of a sine wave in memory and can play back that sine wave as a continuous waveform with very high frequency resolution. But this needn't be a sine wave: you could have any waveform in memory and adjust the rate at which it plays back in the same fashion. That said, I believe these posts will be helpful where I provide more details on how an NCO works:

Numerically Controlled Oscillator (NCO) for phasor implementation?

Problem in understanding DDFS (Direct Digital Frequency Synthesizer.)

As a simplest explanation the approach is to count through the the waveform at a fixed rate (the sampling rate) with an increment that is proportional to the output frequency desired. The higher the precision of your counter (accumulator), the higher precision in frequency can be achieved. If you count with an increment if 1, you will get the lowest non-zero frequency out, and it will take $2^n$ counts to get through the entire waveform in memory once where n is the precision of the accumulator. Only the MSB’s of the accumulator are needed to address the memory based on how big the memory is. If you count by 2 that frequency will twice as fast, etc. So ultimately the frequency resolution is

$$f_\Delta=\frac{f_{clk}}{2^n}$$

Where $f_{clk}$ is the sample rate and n is the number of bits in the accumulator.

And the output frequency assuming the memory contains exactly one cycle of your waveform will be

$$f_{out}= Nf_\Delta$$

Where N is the count increment.

So for example with a 10 MHz clock (update rate) and a 32 bit accumulator which are both typical and practical values, you can generate output frequencies in the 0 to 5MHz first Nyquist zone with 2.33 mHz resolution!  (If combining with an D/A converter the practical upper limit would be approximately 40 MHz to allow for image filtering and an inverse Sinc filter would be considered to compensate for the Sinc reconstruction droop).

The links above go into full details of design parameters including signal quality (SNR and spur levels) and considerations on how many address bits should go to memory etc.


Synthesize Digital Filter with Arbitrary Frequency-Response fdesign.arbmagnphase is a MATLAB built-in function that will create a filter that best fits the frequency-response that it is given.

I'd like to implement a similar function in Python. What algorithms are used to synthesize filters based on an arbitrary frequency response?

For FIR filter design, common algorithms include Remez exchange or Parks-McClellan :
What is the difference between Remez exchange and Parks-McClellan filter design?

For IIR filters, one possibility is the Differential Evolution algorithm, as described in some IEEE journals.  If you have a realistic desired phase response and can guess the filter order, then Berchin's FDLS method may be suitable.


BER vs SNR simulation for GFSK modulation (Closed-form BER for GFSK) I'm in the process of simulating the performance of a GFSK modem in GNU Radio. I have managed to get a curve similar to the curves I see in a number of literature sources. However, I need an accurate baseline that can show how accurate (or far of) the performance is. This might be trivial but I have spent quite a number of hours trying to find out a closed-form BER vs SNR for GFSK. All I have seen so far is BER vs SNR for GMSK only. So my questions are:

Does the closed form of BER for GFSK exist at all?
Is there any explanation as to why it doesn't exist?

All answers and references will be highly appreciated.

Regards,

M.

Well, GFSK is a whole class of modulations, essentially with your freedom being the choice of symbol duration/pulse shaping filter breadth ($\sigma$ in the Gaussian); of this whole class, M makes most sense, so GMSK is really what you'd most probably would want to deal with.

Now, a BER curve really only makes sense if you can describe the "theoretical" receiver of your signal.

In that sense, you could understand GFSK as "extended" GMSK, where you don't use the shortest possible symbol duration, but a multiple thereof, and hence, would need to adjust the $E_b$ in your $E_b/N_0$ curve. But: there should really just be a factor between the $E_b$ in GMSK and generally GFSK, so if you just shift your curve, you'd get a measure of how well your system compares.


How do I zero pad time samples if I only have frequency samples? Related to this question:
Convolution in frequency space with fft and ifft -> num output samples?

When doing convolution in time domain by doing multiplication in frequency domain, you apparently can zero pad the time domain samples if you want the time domain convolution result (ifft of the fft multiplication result), to make the resulting time domain convolution be the "correct" length (len(a)+len(b)-1).

What I'm wondering is what do you do if you don't have the time domain samples yet?

I'm generating two sets of IFFT bin data, doing convolution by multiplying in the frequency domain, and then using IFFT to turn the convolution result into time domain data.

When i do this, the result is the same length as my IFFT bins though, instead of the length i would expect.

How do I lengthen the result in the frequency domain?

Thanks!

assuming I understand it, your question is interesting.  I'm "shooting from the hip" here.  I wonder if it would be useful to you to interpolate your two freq-domain sequences by a factor of two before the multiplication in the freq domain, and then performing your IFFT of the product?
[-Rick-]


Higher-order Kronecker product I am trying to generate a 2D DFT matrix in matlab, which I need for 2D compressed sensing (CS) problems.

Lets say $N_1=8$, $N_2=16$, then according to the requirement of CS, to generate a 2D DFT matrix, we need to first generate DFT matrices  $A$ $(N_1\times N_1)$ and $B$ $(N_2\times N_2)$. Later calculate Kronecker product of $A$ and $B$ which will be of order $N_1\times N_2\times N_1 \times N_2$.

A= dftmtx(N1);
B= dftmtx(N2);dft2D = kron(A,B);

It can be done by Matlab, but there is a limit to the max size one can calculate (in my case its $N_1=64$, $N_2=256$).
But I need it for higher orders i.e., $N_1=128$, $N_2=512$ or higher orders. Is there any efficient way to do this without memory issues?

Since you are using Matlab, I suggest using SPGL1 as your reconstruction solver. It can handle complex numbers and you can specify your matrix operator as a function.

This example demonstrates how to specify your matrix operator as a function instead of a matrix (see opA).

In your case, I would define your operator as a function that:

takes a vector x as input,
reshapes that vector to a matrix X of your desired shape,
applies fft to the matrix: Z = fft(X) (Matlab's fft transforms each column separately if you give it a matrix),
applies fft to the transpose of the result: W = fft(Z'),
reshapes W back into a vector w and returns that

If your operation additionally includes some sub-sampling (which it usually does in order to be compressed sensing) then that can be implemented as indexing into w for the entries that are sampled.

The above steps imply that x is in the domain of the solution you are seeking and therefore sparse. Thus, the sparse vector here is in the time domain and you are sampling in the Fourier domain. If, on the other hand you want to model a signal that is sparse in frequency and you sample in the time domain, you should replace fft by ifft above.


Creating Periodic Impulse Train In MATLAB If we multiply a Continuous-Time signal f(t) by a periodic pulse sequence we can get a sampled version of our signal. The problem I am having is that I am unsure how to create an impulse train in MATLAB that is sampled at fs intervals.

My Attempt at solution (code):

Ts = 1/8000
len = length(inputSignal)
impulseTrain = ones(1,len)
n = 0:Ts:Ts*(len-1)

What I need to do:

conv(inputSignal,impulseTrain)

Okay,so if I was to plot this result it would be very simple however, what I am required to do is convolve this impulse train with an input signal. I am unsure how to pass in the data. How do I get 23000 vectors of length 1 shifted by a distance Ts?

When you load a signal, it is already sampled. So continuous-time does not make much sense here. Creating an impulse train and convolving it cannot help you either. You should know the sampling rate of the signal you are going to load. Assume it is originally sampled at the rate $f_{s,\text{old}}$. This means the samples are apart in time by $1/f_{s,\text{old}}$ seconds.

Now you want to sample this signal with a sampling rate of $f_{s,\text{new}}$. That is, in the sampled signal you will have $1/f_{s,\text{new}}$ seconds time difference between every two consecutive samples. This problem is called resampling, or sampling rate conversion. There are different techniques to do so which you can learn in the given link or similar sources.

A quick approach to do it in Matlab is as follows:

old_fs = 11200; % or any other value it is
old_time_axis = (0:length(input_signal)-1)/old_fs;
new_fs = 8000; % new sampling rate
[p,q] = rat(new_fs/old_fs); %find the ratio of the sampling rates
resampled_signal = resample(input_signal,p,q);
new_time_axis = (0:length(resampled_signal)-1)/new_fs;

Notice that at sharp signal edges distortion can appear. You can minimize such effects, depending on your signal, by adjusting the resampling filter. Documentation can be found here.


Kalman Filter with single input? Is a Kalman Filter still a Kalman filter if there are no inputs that explicitly predict the next state of the system?
For instance, if I have a signal where I know it has characteristics of a first order dynamic system and I have a filter configuration that can in effect discount or ignore changes it considers unlikely is this essentially a Kalman filter without explicit state estimation?

I have tried to follow the derivation for the Kalman filter and I understand that you use some prior knowledge or knowledge of how other system elements may effect your output and so you predict the next state.

Is a sufficient prior to be able to say

  this signal is from a 1st order dynamic system with a time constant of $\tau_x$?

The Kalman filter is the least-squares optimal solution for the states of a set of linear dynamic equations under a Gaussian additive noise assumption. Hence, if you have an unforced linear first-order differential equation you can you use the Kalman filter to estimate the state of this differential equation. I'm not sure it will be very useful though, since the convergence rate of the observer has to be much faster than the dynamics of the first-order system in order to capture most of the transient response, which is the only response in an unforced system,  and to converge if the system is unstable (unless you know the initial state, which is rare). High gain amplifies noise, so you might end up with a poor estimate. For such a simple system you might want to use bandwidth-limited numerical differentiation instead to find the inverse of the filter/system to find the state.


Reading the amplitude values of a FFT in ImageJ Readings the output of the FFT of the picture of a circle in ImageJ at some random point $(163,128)$ I get $X_{\text{Re}}=-0.182$ and $X_{\text{Im}}=-3.854$, using the Complex plot (yellow dot). The value on the output  on the FFT plot is $32$ (cross cursor symbol):

I don't understand why the relation between this latter value ($32$), which I presumed it was the amplitude, and the real and imaginary components, i.e. $\vert X[k]\vert=\sqrt{X_\text{Re}^2+X_\text{Im}^2}$, is not fulfilled.

I had come across the power spectrum logarithmic quote in the ImageJ manual:

The frequency domain image is stored as 32-bit float FHT
attached to the 8-bit image that displays the power spectrum.

but I can't find a formula to translate the values between the output on ImageJ and the square root. For instance $\log 32 =3.466$, while $\sqrt{(-0.182)^2 + (-3.854)}=3.858$. I wonder if there is some normalization across the image, for example to explain this difference.

This is likely related to the Hartley transform used in ImageJ (page 39).

Looking at the source code, the power spectrum image is using a logarithmic scale, adjusted based on the image's dynamic range such that the minimum value (not less than 1) gives a value of 1 and the maximum value gives a value of 254:

if (min<1.0)
    min = 0f;
else
min = (float)Math.log(min);
max = (float)Math.log(max);
scale = (float)(253.0/(max-min));

...
if (r<1f)
    r = 0f;
else
r = (float)Math.log(r);
ps[base+col] = (byte)(((r-min)*scale+0.5)+1);

Making a broad assumption that at least one of those darker pixel is a null would imply that the computed min is close to 0. With this assumption we can proceed to find max as follows:

$$
\begin{align}
   32 &= \left\lfloor \ln(3.858) \frac{253}{\ln(r_\max)} + 1.5 \right\rfloor \\
   31 &\approx \frac{253 \ln(3.858)}{\ln(r_\max)} \\
   \ln(r_\max) &\approx \frac{253\ln(3.858)}{31} \\
   \ln(r_\max) &\approx 11.02 \\
   r_\max &\approx 61020
\end{align}
$$

Since most of the image is white, I would expect a large spike at the 0 frequency point with an amplitude near $NM = 256^2 = 65536$ (where $N$ and $M$ are the width and height of the image). So $r_\max \approx 61020$ is quite plausible.

So for this image (and this image only, since different images would have a different dynamic range), the approximate scaling would be:

$$
\begin{align}
\mbox{value} &\approx \left\lfloor 22.96 \ln(\sqrt{X_\text{Re}^2 + X_\text{Im}^2}) + 1.5 \right\rfloor \\
  &\approx \left\lfloor 11.48 \ln(X_\text{Re}^2 + X_\text{Im}^2) + 1.5 \right\rfloor
\end{align}
$$


Sampling the impulse response Given two IIR filters (one Low-Pass and one High-Pass) designed using Butterworth, how can I know if they can be obtained using sampling the impulse response?

The short answer is: You can implement a lowpass filter by sampling its prototype continuous time impulse response $h_c(t)$, but you cannot implement a highpass filter by sampling its prototype impulse response, due to highpass filter being  non-bandlimited i.e. there will be unavoidable aliasing due to the definition of the highpass filter. But read on.

The technique of obtaning a discrete time filter impulse response $h_d[n]$ from a continuous time prototype filter by sampling its impulse response $h_c(t)$ as given by $$h_d[n] = T_s h_c(nT_s)$$ is called as the method of impulse-invariance , where $T_s$ is the sampling period.

This operation of obtaning a discrete time filter from a continuous time prototype is analogous to the more traditional operation of signal sampling. In fact their implications are the same: There will be aliasing at the resulting discerete time spectrum of the converted filter if either 1-The continuous time filter is not bandlimited or 2-The sampling rate is not adequate, given by the Nyquist rate as applied to sampling of bandlimited signals.

A continous time lowpass filter by definition, is always bandlimited, and we can always, therefore, find a small enough $T_d$ which will avoid aliasing at the resulting discerete time spectrum of the filter being converted.

On the other hand, a continuous time highpass filter is going to a have an infinite bandwidth and as such the sampling period required is $T_s = 0$ , i.e. it cannot be sampled.

However, if you will be carefull enough, and if you can be sure about the bandwidth requirements of every stage in your signal processing chain, then may be you can replace a highpass filter by a bandpass one, when such is the case, then you can still apply the technique of impulse invariance to obtain a discrete time band pass filter from a continuous time bandpass filter prototype, by choosing a suitable sampling period $T_s$

Note that, the method of impulse invariance produces an exact linear mapping from continuous time to discrete time frequency spectrum. It will preserve all the characteristics of the continuous time filter (unlike the method of bilinear transform which warps the spectrum)


How to match two signals that have similar shape but are scaled and expanded? I have two signals A and B in time domain as follows:

Visually both signals look very similar; however the X and Y axis range are different for both the signals. Signal B is basically a "stretched" version of signal A. The minimas, maximas occur at different x and y values; but relative distance between successive peaks of the signals have somewhat an uniform proportion.
I'm looking for a similarity metric that would give high similarity value for such pairs. Note dimensions of both signals are different. Due to the domain constraint, I cannot pad values to any of them. I tried cross correlation, circular correlation, and other common similarity measure but they didn't give satisfactory results.

You can always match the extrema via a simple linear mapping, e.g. $y[n] = \frac{x[n]-min(x[n])}{range(x[n])}$ where $range(\cdot)$ is the absolute difference between the maximum and minimum values and $x$ is the input signal. This will get $y \in [0 .. 1]$ and then you can shift that interval to whatever (min,max) and baseline you like.

The dilations in time might be a bit harder to deal with.

If you are positive that the dilations are more or less by a constant factor, then prior to assessing the similarity between the two waveforms you can resample the longer one down to the length of the shorter one and then apply correlation. In fact, you could still base your similarity only on the comparable portions of the signals spectra, either by a direct error term or via something like coherence.

But if the dilations are not stable, then Dynamic Time Warping (DTW) (as was suggested) is the better option to go for as it will give you a "speed up / slow down" vector which you can use to either map one signal on top of the other or directly, as a similarity / disimilarity metric.

Hope this helps.


Characterization of the amount of modulation in a signal We have used the ICRA babble noise to generate an 8-talker babble. As a background - babble noise is used as a substitute for competing talkers as babble noise lacks lexical content, but contains the same envelope and frequency content as 'real' speech. Babble noise is, therefore, less distracting for listeners and hence it's a more direct, 'peripheral' type of masker. Lexical masking occurs more centrally.

Anyway - we are interested in whether the masking effects we see are due to the slowly varying part of the babble noise (the envelope) or the frequency content.

Now, because we used an 8-talker babble using asynchronous, random starting points, the amount of modulation in the masker is less, as all the 8 babble noises are randomly added. Hence, it becomes more of a static type of noise, and  listening in the gaps is hardly possible. The latter is especially important - are there gaps present in the noise, and if yes, how long are they? Likely there are no true gaps anymore, but only stretches of noise where relatively little noise is present, so shallow gaps in the noise.

Now, I wish to compare the characteristics of the single-talker original noise and our 8-talker babble in terms of:

amount (number/density) of gaps present in the noise
the gap duration
and amplitude of the gaps (modulation depth)

I think a good starting point is to determine the envelope of the fluctuating noise using a Hilbert transform. But then, how do I go about analyzing the envelope in terms of the two parameters, namely gap-duration and amplitude? I have added the envelope of a shorty stretch of the original single-talker babble below (x-axis label is seconds):

FWIW:

I'm used to do my signal analysis in Matlab;
Noise files were software-corrected for rms for overall sound level;
Noise was played back in a sound-proof booth and we will record the signals back with a KEMAR manikin to include the head shadow effect.

As far as the equivalent of an "envelope detector", consider using a moving average of the absolute value (as it would be simple to implement, if you like the approach more elaborate filters could be used). This is how you would do that in Matlab:

M = 10
out = filter(ones(M,1), 1, waveform)

Where M is the number of samples that you want to actually do the moving average over. Note that you will be making a distinction between what you consider the frequency components of your signal vs the frequency components of the envelope as this algorithm seeks to distinguish the two, with the envelope being a low frequency component.

To determine gap duration, you can simply set a threshold on the envelope based on what you decide to be a gap or not and then create a bipolar two level signal (+1 or -1) to measure the duration of the gaps:

gapTest = abs(out- threshold)

As for setting the threshold itself, I would be tempted to use an "N-Sigma" appraoch, where the threshold is determined based on the standard deviation of the data (so would scale automatically and be consistent across sets).

And the following will report the length as "durations" and quantity of gaps as "count" (note that I did it such that leading and trailing gaps are ignored but you may decide to do that differently):

# find edges
x = diff(gapTest)

# ignore leading gaps (as they are partial)
fall = find(x == -1)
x2 = x(fall(1):end)

#find gaps and duration
rise = find(x2 == 1)
fall = find(x2 == -1)

#ignore trailing gaps (as they are partial)
len = min(length(rise), length(fall))
durations = rise(1:len) - fall(1:len)
count = length(durations)

You will need to test the above for all cases and modify accordingly.

As far as the equivalent of a "mod depth" consider using the mean of the data above your threshold compared to the mean of the data below your threshold.


Filter for Reducing Background Noise in Image I have images where there is a lot of "black" background (few shades of black). (In many images at least half the pixels are background). I need to get interest points from the image, but because the background is not smooth (although black) , I get a lot of "false" key points in the background.
What filter/s will fit best to keep the details in the image but will smooth the background as much as possible ?
I tried median filter, but the size needed to reduce background noise ([5 5]) was to strong for the small details in the image I wanted to keep sharp, same for Gaussian with very small sigma.

Thanks

You can reduce image size or use scale-invariant interest point detector which can deal with blurry images.

Please note that using some elaborate denoising algorithm may change each image in a slightly different way, thus possibly distorting the positions of your interest points.


doubt of intersection of ROC Z Transform dear friends of StackExchange.
I have a doubt of the intersection of two ROC.
I have H(Z), X(Z) and and i have to determine:

$$ \begin{align} Y(Z)= H(Z)X(Z)\end{align}$$
$$
\displaystyle
$$

$ ROC $ H(Z)  $ \cap $ $ ROC $ X(Z)
$$
\displaystyle
$$
poles of $H(Z): p1=  \left |  -\frac{1}{4}  \right | $ ; $p2= \left |  \frac{3}{2}  \right |$
$$
\displaystyle
$$

The system is supposed stable. Ops, sorry i forgot to write it in the initial post, under this assumption
the corresponding $ROC$ of $H(Z)$ is :
$$\begin{align} \quad&\ \left | \frac{1}{4} \right | <z<  \left | \frac{3}{2} \right | \end{align}$$
$$
\displaystyle
$$
  pole of $X(Z): p3=  \left |  2 \right |  $ $$ \displaystyle
 $$
For stability condition:
 $$\begin{align} \quad&\ z<  \left |2 \right | \end{align}$$
The intersection area is:
$$ \displaystyle$$

$Y(Z):  ROC $ H(Z)  $ \cap $ $ ROC $ Y(Z) = $$\begin{align} \quad&\ \left | \frac{1}{4} \right | < \left |  z \right | < \left | \frac{3}{2} \right | \end{align}$$

i see it drawing the circumference and tracing the circumferences related to the poles and their corresponding ROC
 .

In order to answer this question we would need to know the ROCs of $H(z)$ and $Y(z)$. The pole radii only determine the limits of the possible ROCs. E.g., for $H(z)$ there are three different possible ROCs:

$$\begin{align}1.\quad&|z|>\frac{3}{2}\\
2.\quad&\frac32>|z|>\frac14\\
3.\quad&|z|<\frac14\end{align}\tag{1}$$

For $Y(z)$, there are two possible ROCs:

$$\begin{align}1.\quad&|z|<2\\
2.\quad&|z|>2\end{align}\tag{2}$$

Note that each of these ROCs corresponds to a different time domain sequence.

Certain combinations of the ROCs in $(1)$ and $(2)$ overlap, others don't. Below are shown the combinations that result in an overlap of ROCs, and the resulting combined ROC. The first number is the number of the ROC of $H(z)$ given in $(1)$, the second number is the ROC of $Y(z)$ given in $(2)$:

$$\begin{align}1+2:\quad&|z|>2\\
1+1:\quad&2>|z|> \frac32\\
2+1:\quad&\frac32>|z|>\frac14\\
3+1:\quad&|z|<\frac14\end{align}\tag{3}$$

All other combinations of ROCs result in no overlap.


why 2 vectors are not the same in matlab? I have 2 vectors:

a = 1.0e-04*[0.110000000000000   0.140000000000000   0.140000000000000   0.140000000000000];
b =         [0.000011000000000   0.000014000000000   0.000014000000000   0.000014000000000];

by look at by eyes, it seems similar; but when I use this command:

a == b

then I got:

ans =
  1×4 logical array
   0   0   0   0

The expectation is:

ans =
      1×4 logical array
       1   1   1   1

Does anyone knows what wrong with this command in the matlab?

Thank you.

Neither 0.11 nor 1E-4 can be exactly represented in floating point arithmetic. And the product of 2 numbers will introduce an extra error.

Testing for equality in floating point arithmetic is therefore more complex than it seems.


Selecting gradient orientation for Generalized Hough Transform while scaling and rotation is present I am using this link to study Generalized Hough Transform. I will explain it briefly. In order to detect a shape, an R-table is created for the shape with $\phi$ indicating gradient orientation for points on the edge of the shape. Then $(r, \beta)$ is used to indicate the vector from each of these points to the centroid in the shape. While detecting the shape, we will first find all the edge points and their gradient orientation $(\phi)$. Based on this $\phi$, $(r, \beta)$ values are used to calculate centroid points. Based on the votes on centroid points (and R table), the shape can be detected.

My question is related to rotation of the shape which is explained towards the end. If the shape is rotated by $\theta$, then the gradient orientation $(\phi)$
 for a given edge point changes. So, shouldn't we do either one of the following:

1) Rotate all the $\phi$ values in R table by $\theta$?

OR

2) Rotate gradient vector by $(-\theta)$ and then calculate $\phi$ for the edge point?

If the shape is rotated by $\theta$, then the gradient orientation ($\phi$) for a given edge point changes. So, shouldn't we do either one of the following:

  Rotate all the ϕ values in R table by θ?

  OR

  Rotate gradient vector by ($-\theta$) and then calculate $\phi$ for the edge point?

The generalised Hough transform does #1 implicitly by the $S, \theta$ parameters when accumulating.

Note here that usually, you don't really know that the extracted contour has been rotated by some angle (let's call it $\alpha$). You usually are trying to find that angle.

During the preparation stage, you construct the "template", $\phi$. This quantises the circumference of a circle that encloses the shape to be recognised into $k$ levels ($k$ distinct angles). Then, for each $k^{th}$ level, you look at the contours that make up your template image at the $\phi_k$ bearing and note down their $(r, \beta)$. The $(r,\beta)$ is nothing more than the polar coordinates of an element of your template image's contour that is oriented at angle $\phi_k$.

If you wanted to recognise a simple circle, then you would only have one column of $(r, \beta)$. If that circle enclosed another, smaller circle (imagine for example a road sign), then at some orientation $\phi_k$ you would have more than one "edges" and therefore more than one columns of $(r, \beta)$ data.

Now, during the recognition phase, you use the $\angle G$ to lookup an entry in $\phi$. Once you have found that entry, your take that entry's $(r, \beta)$ and use them to find the "right Hough counter cell" $H(x_c,y_c)$ and increment it.

The point here is that the $(r, \beta)$ are calculated with respect to the centroid and it is that centroid $x_c,y_c$ that accumulates quickly in $H$. When the shape is matching, it means that your test image contains a lot of $\angle G$ entries that have resulted in the $H(x_c,y_c)$ counter getting a very high value.

The next step is to be able to recognise the template $\phi$ irrespectively of scale $S$ and orientation $\theta$.

Notice here that if you were to take the Hough transform of the template with its centroid translated at $(0,0)$, then different rotations $\theta_1, \theta_2, \theta_3, \ldots$ are basically cyclic lateral shifts in the Hough space (where the $x$ dimension now corresponds to bearing $\theta$) and different (uniform) scalings $S_1, S_2, S_3, \ldots$ correspond to different "arc amplitudes" (in the $y$ direction that corresponds to distance from origin).

So, the template doesn't really change, all that you do is consider different scalings and rotations at the same time. From another point of view it is as if you were running $N$ different Hough transforms where you scale / rotate the original image by $(S,\theta)_n$. You would produce one $H$ for every $(S, \theta)_n$ such as $\left(1,0\right), \left(1, \frac{\pi}{4}\right), \left(1, \frac{\pi}{2}\right), \left(1, \frac{3 \pi}{4}\right), \left(1, \pi\right), \left(0.5,0\right), \left(0.5, \frac{\pi}{4} \right), \ldots, \left(0.25, 0 \right), \ldots$. In this case here we are looking for 4 orientations at scale 1, 0.5, 0.25 and so on.

But, rotating and scaling the original image can be costly because there are definitely more pixels in an image than elements in the $k$ level quantised periphery of the curve you are trying to find.

For this reason, instead of rotating the original image and running parallel Hough transforms, you simply scale and rotate $\phi$ but you do it on the fly. Instead of modifying the $\phi$ you just introduce $S, \theta$. $S$ modulates your $r$ (the radius of the template) and $\theta$ offsets your $\beta$ (the bearing of your template).

In this way, you produce $N$ matrices $H(x,y)$, one for each $(S, \theta)_n$ you would like to consider. Again, the $x_c, y_c$ counter that crosses a certain threshold marks the location of the recognised contour but in addition now you can also infer at which scale and rotation by fixing the $(S, \theta)$ pair.

Hope this helps.


How does 3D Time-of-Flight (TOF) camera work? This is my first post in this forum! So Point me if anything is not clear!

I want to know how TOF camera works? So I've done Google & read some nice explanation about How Phase difference between Transmitted & Received signals can be used to calculate the Depth of an Object from Camera with good example.

What i did understand:

Depth Can be measured by Calculating the Phase change of the Transmitted signal.
How a light source(IR) can be modulated & demodulated ? Answer
How to filter the reflected IR signal having the same modulation frequency? Answer

What i cannot understand:

How phase change of a continuous square wave signal can be used to Measure the depth?

Update: Here is a Nice Explanation about this.

Note: I'm talking about TI & Mesa Kind of TOF Camera's using IR Led

The (unwrapped) phase between two constant same frequency signals is linearly related to any delay between the two signals.

Measuring phase of a fixed frequency periodic signal over many periods is sometimes more accurate than trying to detect delay locally (by zero crossing, etc.)


How do you change pitch with wavelet transforms in MATLAB? I'm trying to increase the pitch of a wav file using MATLAB. (I don't expect it to sound very good, but I'm trying to figure out how wavelets in MATLAB work.) Here's what I'm doing:

scales = {2.34375,0.5,21};
result = cwtft(b5,'scales',scales);
result.scales = 2*scales;
b5high = icwtft(result);

(b5 is just the opening of Beethoven's Fifth, but it could be any sound file. It happens to have a sample rate of 48 kHz.)

To my thinking, I'm doing a wavelet transform of b5 using a Morlet wavelet. The scales that I'm using capture the human-audible frequencies for a 48-kHz signal. I then double the scales vector--resulting in a shift upward of 1 octave, and then convert back into a signal. However, when I play b5high, it sounds exactly like b5 - no pitch shifting at all.

What's wrong with my logic?

Thanks.

I haven't implemented a pitch shifter using wavelets, but I am currently learning about them. From what I understand scales are the inverse of frequency, therefore to get your desired result you'd have to divide the scales instead of multiplying.

Furthermore, according to 'Applications of the Continuous Wavelet Transform in the Processing of Musical Signals', you also have to update the phase as adjusting the scales isn't enough. The algorithm outlined in the paper is:

semitone = 12;
pitchFactor = 2^(semitone/12);

wt = cwtft(input);

mag = abs(wt.cfs);
phase = angle(wt.cfs);
unwrapped = unwrap(phase);

wt.cfs = mag.*exp(1i*unwrapped*pitchFactor);
wt.scales = wt.scales./pitchFactor;

result = icwtft(wt);

Sadly, that's all I can offer until I learn more.


Calculating the complex signal's average power $$x(t) = \cos(\pi i t/20+\pi/4) - 2je^{j \,12\pi i t} + 5\sin(2\pi i t/3+\pi/3) $$

I need to find the signal's average power.
How can I do that?
The $i$ index (represented in red in the below image) just denotes a random integer.

Solve the integral:
$$
P = \lim_{T \to \infty} \frac{1}{T}\int_{-T/2}^{T/2} |x(t)|^2 dt
$$
This is usually unwieldy to solve directly; it can be shown that if a signal is periodic, you only need to integrate over the fundamental period $T_0$:
$$
P = \frac{1}{T_0}\int_{-T_0/2}^{T_0/2} |x(t)|^2 dt
$$

EDIT:

Another "trick" is to note that sinusoidal signals of different frequencies are incoherent, so when integrating over long periods of time, it is only necessary to find the power of the individual components because the cross terms go to zero.

This means that
$$
P = P_1 + P_2 + P_3
$$
assuming $x_1$, $x_2$, and $x_3$ are incoherent.  Find the individual powers and sum them.


What do you call the random Gaussian vectors in compressed sensing? Let $(e_i)_{1\leq i\leq M}$ be vectors with zero-mean i.i.d. entries from a Gaussian distribution. In a compressed-sensing setup, I have observed the collection of scalar products $(\langle x\, |\, e_i\rangle)_{1\leq i\leq M}$ for an unknown vector $x$.

In the compressed sensing literature, what is the name for the collection of vectors   $(e_i)_{1\leq i\leq M}$?

The collection of scalar products is often called measurement. The collection of vectors (each being an atom) can be called measurement dictionaries or sensing dictionaries. This is not specific to random atoms.


A problem about whitening of a signal I intend to whiten the spectrum of a time series using a point-wise normalization in the frequency domain.

In my first test, it looks pretty good.

figure;
n = 2e3;
t = 1:n;
m = n/2+1;
f = linspace(0,1,m);
x = rand(1,n) - 0.5;
y = fft(x);
subplot(2,2,1); plot(t,x,'k');
axis tight; xlabel('t/s'); title('original signal x');
subplot(2,2,2); plot(f,abs(y(1:m)),'k');
axis tight; xlabel('f/hz'); title('original spectrum');
y1 = y ./ abs(y);
x1 = real(ifft(y1));
y2 = fft(x1);
x2 = real(ifft(y2));
subplot(2,2,3); plot(t,x1,'k',t,x2-x1,'r');
axis tight; xlabel('t/s'); title('whitened signal x1(black) and x2-x1(red)');
subplot(2,2,4); plot(f,abs(y1(1:m)),'k',f,abs(y2(1:m)),'r');
axis tight; xlabel('f/hz'); title('whitened spectrum y1(black) and y2(red)');

The spectrum is flat after whitening.

But when I windowed the spectrum, sth strange happened.

figure;
n = 2e3;
t = 1:n;
m = n/2+1;
f = linspace(0,1,m);
x = rand(1,n) - 0.5;
y = fft(x);
subplot(2,2,1); plot(t,x,'k');
axis tight; xlabel('t/s'); title('original signal x');
subplot(2,2,2); plot(f,abs(y(1:m)),'k');
axis tight; xlabel('f/hz'); title('original spectrum');
y1 = y ./ abs(y);

k = n/10;
w = sin(linspace(0,pi/2,k));
y1(1:k) = y1(1:k) .* w;
y1(m:-1:m-k+1) = y1(m:-1:m-k+1) .* w;
y1(n-m+3:n) = y1(m-1:-1:2);

x1 = real(ifft(y1));
y2 = fft(x1);
x2 = real(ifft(y2));
subplot(2,2,3); plot(t,x1,'k',t,x2-x1,'r');
axis tight; xlabel('t/s'); title('whitened signal x1(black) and x2-x1(red)');
subplot(2,2,4); plot(f,abs(y1(1:m)),'k',f,abs(y2(1:m)),'r');
axis tight; xlabel('f/hz'); title('whitened spectrum y1(black) and y2(red)');

As shown in subplot 3, x1 and x2 are the same (the difference is zero; see red line), but their spectrum in subplot 4 (black and red) is different.

After IFFT y1 to x1 and FFT x1 to y2, y2 differs from y1.

As suggested by @Drazick, I post as an answer.

The good code according to the comments:

figure;
n = 2e3;
t = 1:n;
m = n/2+1;
f = linspace(0,1,m);
x = rand(1,n) - 0.5;
y = fft(x);
subplot(2,2,1); plot(t,x,'k');
axis tight; xlabel('t/s'); title('original signal x');
subplot(2,2,2); plot(f,abs(y(1:m)),'k');
axis tight; xlabel('f/hz'); title('original spectrum');
y1 = y ./ abs(y);
x1 = real(ifft(y1));
y2 = fft(x1);
x2 = real(ifft(y2));
subplot(2,2,3); plot(t,x1,'k',t,x2-x1,'r');
axis tight; xlabel('t/s'); title('whitened signal x1(black) and x2-x1(red)');
subplot(2,2,4); plot(f,abs(y1(1:m)),'k',f,abs(y2(1:m)),'r');
axis tight; xlabel('f/hz'); title('whitened spectrum y1(black) and y2(red)');

The spectrum is flat after whitening.


How would sigma-delta modulation work in software? I've recently come across this YouTube video that shows a PCM stream being played back via the 1-bit internal "PC Speaker", and after hours of googling I've come to think this is done via Sigma-Delta modulation.

So far, what I gather from the process is that it takes the original digital signal, oversamples and approximates it via PWM while applying a low-pass filter to the 1-bit output. I'm not quite sure if this is right.

Not having a background in DSP, I'm wondering:
How could this process be implemented in software?

Upsample the signal to desired frequency. You can use linear interpolation (poor quality but simple and fast) or sinc interpolation (good quality but slow). For example of library dedicated to good quality resampling, see libsoxr: http://sourceforge.net/p/soxr/wiki/Home/
Perform dithering with noise shaping. You need to filter the dithering signal to minimize noise within 0-20000 Hz range while increasing it in ultrasonic range. For details on noise shaping, see this excellent document: http://www.beis.de/Elektronik/DeltaSigma/DeltaSigma.html especially Figure 9 which shows adding noise in feedback loop.

That's how delta-sigma DACs work and this can certainly be implemented in software.

PWM would require upsampling 65536× for 16-bit audio so it isn't used. PDM (Pulse Density Modulation) is used instead and is essentially the same as delta-sigma modulation. And it requires only 64× upsampling for Hi-Fi audio.


Frequency Spectrum of a sequence The OQPSK modulator in Matlab produces an output sequence as follows:
[0.707
0.707+0.707i
-0.707-0.707i
...and so on
]
This is the sequence and I am trying to view the frequency spectrum of this sequence. With the sampling rate being 2 MHz, I expect to see a peak centered at 0 with a bandwidth of 2 MHz. When I upsample the above sequence and multiply this sequence with a carrier of say 2.4 GHz, I should have 2 peaks centered at -2.4 GHz and 2.4 GHz and zeros at other frequencies. But somehow, the spectrum is distributed between -2.4 and 2.4 GHz. Am I making a mistake in plotting the frequency spectrum or is there something wrong in what I understood?

You're missing a step: pulse shaping. The sequence [.7+.7i, ...] is a sequence of pulse amplitudes. As such, it should not be considered to be signal; in consequence, it does not have a spectrum.

Say your sequence of amplitudes is $A=a_0, a_1, a_2,\ldots$. Now, you need to specify your symbol rate (how many pulses per second you'll transmit). Call this rate $R_p$ and let $T_p=1/R_p$ be the pulse interval. Finally, you need to specify a pulse $p(t)$ that has this property: $$\int_{-\infty}^\infty p(t-kT_p)p(t-lT_p)dt=\begin{cases}1,\textrm{ if $k=l$}\\0,\textrm{ if $k\neq l$}\end{cases}$$for $k$ and $l$ integers. Now you're ready to create your actual signal: $$s(t)=\sum_i a_ip(t-iT_p).$$ This signal:

Has a spectrum; specifically, its bandwidth is the same as the bandwidth of the pulse $p(t)$.
Can be up-converted to a carrier frequency.
Can be converted to an analog signal and physically transmitted over a medium.

For completeness, I'll explain how to recover the amplitudes $a_i$ from $s(t)$ (ignoring noise for simplicity). The receiver can estimate $a_k$ by doing this operation: $$\hat{a_k}=\int_{-\infty}^\infty s(t)p(t-kTp).$$ You can easily verifty that this is the case, using the properties of $p(t)$. Such an operation can be implemented with a matched filter.


What kind of modulation is this? I'm trying to find out what kind of modulation and encoding is used by my window cover remote control, aiming to implement the same on a Raspberry Pi to add a bit of home automation.

A colleague helped me sample the signal with an USRP and gave me the data files. I have examined them with Inspectrum.

I'm afraid the samples have been recorded centered very close to the frequency of the signal. I've later read that this might not be a good idea, and I think this is what causes the permanent noise you can see in the pictures. Should I center the sampling a few kHz below in the future?

Single press of the "open window" signal

Symbol rate, or rather what I guessed was the symbol rate, seems consistent at this level (around 42 bauds).

This shows frequency energy. X axis is time, Y axis is frequency (baseband). Colors have been tuned to show only the interesting power.

Detail of one of the "blocks"

Here, the last picture shows, frop to bottom: frequency, IQ (red/blue), and phase. I don't quite know what the frequency chart represents here.

Zoom on the transition between sparse to dense

The "symbol rate" doesn't seem constant here.

I have no idea what kind of modulation, encoding, bit rate, or protocol this remote uses, and right now my only option would be to blindly try all the demod blocks in gnuradio.

It'd be extra helpful if the answer could include some guidance as to how to demodulate this signal in gnuradio in order to get bits out. I think I can take it from there.

The controller is for a Velux window cover, but I haven't found any information on the web. I ultimately need to find out what chip to buy for the Arduino/rPi.

I can tell the chip has some text inscribed, but I cannot read it, seems it's worn off, but I don't have a magnifier... :/

Raw data uploaded to: https://www.dropbox.com/s/rh2k7ho68dvoxhd/data_mando3.dat?dl=0 . Sample rate is 3MHz. Data format is GNU Radio Companion default, which is IQ, each component expressed as a 32 bit float.

Update

After some more investigation this seems to be io-homecontrol protocol. It seems it hasn't been reverse engineered, it's ciphered, and no information is available. It's a two way protocol so I'll need to grab samples from full conversations (the current samples are taken with only the remote available).

This chip in theory is able to work with this protocol: http://www.analog.com/media/en/technical-documentation/data-sheets/ADF7022_2page.pdf

The consortium will not provide specs. I'm still interested in sniffing a few conversations with GNURadio and see if I can work from them. Given this, I can only recommend against io-homecontrol and Velux products.

Thank you!

Ok I did some signal forensics on the data capture and believe the modulation is a form of FSK.

The FSK modulation was +/- 20 KHz with a data rate of 38 KHz.

UPDATE:  The OP discovery that this is "io-homecontrol" and the datasheet from ADI that he found has confirmed that this is indeed FSK with a deviation of 20KHz and 38.4 Kbps data rate.  Further this link gives more details on the modulation format:

https://cansecwest.com/slides/2015/From_Baseband_to_bitstream_Andy_Davis.pdf

It looks like it may actually be Gaussian-FSK, with full-response signalling (BT=1), meaning the response for each bit is completed before the next bit starts. This can be easily implemented with a Gaussian filter (see Gaussian FIR filter with no multipliers?) followed by an NCO. In this case each data bit is represented as impulses followed by zeros for the full length of the filter and the output of the filter drives the frequency control word into the NCO (a digital VCO). (For partial-response signalling which requires a more complicated receiver and not what they are doing here, the zeros are shorter than the length of the filter, but that would also work).  Frequency deviation can be set with a multiplier gain constant between the filter output and NCO input. The choice to go with Gaussian FSK vs simple FSK is for better spectral containment in the transmitted waveform (and is the reason for the rounded transitions that we see).

I will include the plots in case anyone has further insight:

First a macro view in vicinity of the first burst, occurring around 0.93 seconds. To get this plot, I had removed a frequency offset corresponding to 0.2688 radians/sample after decimating by 10 samples, assuming the original sampling rate was 3 MHz, this corresponds to a 12.834KHz offset.

This plot shows the unwrapped phase vs time (vertical axis is radians) and the corresponding magnitude vs time. I assume when the magnitude is low that we are viewing the noise, and the signal of interest is occurring only when the magnitude is higher.

Zooming in on the first part of the burst, the pulse starts at approx 0.9265 seconds, and at the start of the burst, the phase is "flat" for approximately 3 ms.  Corresponding to an unmodulated carrier, I believe, from comparing to the detail we see later.

After the 3 ms of the "flat phase" portion, there is a very long pattern of what appears to be 1 0 1 0 1 0 with a data rate of 37.5 KS/s.  The modulation is going between 0 and 180 degrees, but later we will see clear phase ramps making me suspect FSK instead of BPSK. (If it was MSK I would have expected to see evidence of 90 degree rotations, but the minimum rotation I found was 180 degrees).
Note at the very start of the pattern (assuming +180 =1 and -180 =0), we see 1 1 0 1 0 1 0 1... which then lasts with the 1 0 1 0 pattern for approx 90 ms before again inserting a "1 1" near the 1.02 second mark, and then continuing again as 0 1 0 1 0 1...  until the 1.414 second mark. (approx 0.4 seconds duration with a 1 0 1 0 pattern).  Note the "hump" in the overall trajectory; I believe that we are just seeing the wandering of the clock (either the radio or the sampling clock that took the data) given the long time durations involved.

At 1.414 seconds the real data starts. (The first was probably a sync/acquisition pattern).  Looking first at a macro view we see a pattern that at this view appears to repeat 3 times before the burst ends.

Zooming in on the first region, we see consistent positive and negative phase sloped indicative of FSK. At the end of the data burst there is a long duration of 1's for 4 ms followed by the 101010 pattern again for about 14 ms.

The next burst starting around 1.51 seconds, shortly after the first one again has the same starting pattern of an unmodulated carrier followed by the 101010 long modulation pattern.  (Although slightly different in that the start was 1 1 1 1 0 1 0 1 0 1 0 1 but the duration of the 1010 pattern before the data modulation was again 14 ms)


Multidimensional FFT with mixed directions Let's say I have three-dimensional complex array $X[m,n,p]$, $m \in \{0, 1, \ldots, M-1\}$, $n \in \{0, 1, \ldots, N-1\}$, $p \in \{0, 1, \ldots, P-1\}$; i.e. $X \in \mathbb C^{M \times N \times P}$.   I would like to take a three-dimensional DFT of this array; along dimensions $m$ and $n$ I would like to use the forward transform, and along dimension $p$ I would like to use the inverse transform.

Every multidimensional FFT implementation I've looked at (FFTW, MKL, numpy.fft) just has one parameter to specify as "forward" or "backwards/inverse", implying that the same transform is applied in each dimension.

Is there an "easy" way to do this mixed-direction FFT directly without having to break it up into two steps (i.e. doing a 2D forward transform and then a set of 1D inverse transforms)?

Either a conceptual answer (e.g. "conjugate this then do this") or a "programming" answer (e.g. "use this flag in MKL") would be fine.

if I understood it right, given your 3D volumetric data $x[m,n,p]$, first you want to take a 2D forward DFT (via 2D FFT) in the first two variables $m,n$ (i.e., for each m-n plane along 3rd dimension) and then take a 1D inverse DFT of vectors along the third dimension.

Note that in the first operation you replace each data plane with frequency plane and create a new volumetric data, then you will be taking inverse DFT of this resulting volume, the MATLAB/OCTAVE code below demonstrates a simple example:

A = zeros(2,2,3);      % create a 3D data array
A(:,:,1) = [1,2;3,4];  % first plane of 3D data
A(:,:,2) = [1,2;3,4];  % second plane of 3D data
A(:,:,3) = [1,2;3,4];  % third plane of 3D data

B = fft2(A)        % it'll take 2D-DFT independently for each plane of A

C = ifft(B,[],3)  % we indicate that ifft() be taken along 3rd dimension...

As you can see, since after the first step we have same frequency data for each plane of B, the third step produced impulses for each vector of C along 3rd dimension, due to inverse DFT.

I believe this is easy enough? But still in two steps (you can merge them into one line but still two calls)...


Harmonics below oscillation frequency in the brain Apologies if this question is trivial, I'm a neuroscientist and would hugely appreciate your help to correctly interpret my analysis.

I am recording EEG signals and am manipulating the brain using light flickering at 10Hz. I then compute the PSD of this signal using FFT of 4s long epochs. I am interested in what the manipulation does to frequencies around 5Hz.

My problem is that my 10Hz stimulation creates an artifact at 10Hz (light induces current flow). I see this clearly in the PSD and I also see higher harmonics at 20 and 30 Hz.  My question is whether I should expect this artifact to affect the 5 Hz frequency as well? I don't see a clear peak but I don't fully understand the theory so I'm not sure.

You explain in the comments that the signal is sampled at 50 KHz, then low pass filtered at 128 Hz with a 2nd order biquad filter and then resampled and stored at 256 Hz.

If you weren't prudent with choice of sampling rate and low pass filtering, then the harmonics of 10 Hz could create a component at 5 Hz through aliasing. For example this would occur if you used a sampling rate of 25 Hz: the third harmonic at 30 Hz would fold to the 5 Hz frequency. I explain aliasing in more detail at this post here:

Where should I set my anti-aliasing filter corner frequency for this signal?

That said, there are two prominent mechanisms that could create folding images in vicinity of 5 Hz due to aliasing during this sampling process.

The first would be due to interference components above the Nyquist frequency of 25 KHz prior to the sampling of the analog signal - specifically any frequency components at $N f_s \pm 5$ Hz would fold to 5 Hz in your digital signal where $N$ is any integer and $f_s$ here is 50 KHz. (So 50 KHz +/- 5 Hz, 100 KHz +/-5 Hz are example locations where interference can alias in). To avoid this be sure to use a good anti-aliasing filter and pay particular attention to its rejection around $N f_s \pm 5$ Hz with a rejection band around those frequencies that are within the frequency range that you want to maintain to be free of interference for your measurement.  Given this is many of orders larger than your 10 Hz interference, it is less likely this is an issue however should not be discounted depending on your measurement sensitivity desired.

The second mechanism, and likely dominant one in this case, is your resampling approach. Resampling to 256 Hz is a fractional rate conversion from 50 KHz and there could be many opportunities to create harmonic folding effects. The actual implementation needs to be reviewed in detail including significantly the performance of the resampling filter prior to decimation. Resampling to 256 Hz is identical to sampling at 256 Hz directly, so any frequency content above the Nyquist frequency of 128 Hz would fold into band. Specifically in this case the frequencies at $N 256 \pm 5$ Hz would land at 5 Hz in the digital domain.

Without seeing the actual implementation or ability to change it this would be difficult to predict although certainly can be measured. Induce a 10 Hz test signal with high harmonic content (impulses) and measure the results without your actual device you are measuring connected to assess your measurement noise floor in that condition.

Regardless of aliasing, to be sure that the interference itself does not contain 5 Hz components, I would also recommend analyzing a measurement without your true data source connected (such that you are only evaluating background noise). This would also give you a fidelity of your measurement system in showing you what sensitivity can be achieved, and the related confidence of your results.


What does it mean to rotate an image by 0.1 degrees? In the old Tangent Prop paper (you don't need to know this paper), the tangent vector of the rotation transformation (think of this as a vector-valued function) at some input image is approximated by (1) rotating this input image by a very small amount, (2) subtracting this rotated image from the original image and (3) dividing the result of subtraction by the degree rotated.

However, I'm not understanding what it means to rotate an image (i.e., a matrix) by a very small amount (less than one degree).

Also, how can this be implemented in Python using libraries like Python / Scipy / Skimage?

In addition, how am I supposed to know that the image has actually been rotated, since the rotation is so small?

Note: I tried to do rotation with Skimage, but it seems like Skimage allows rotations of at least 1 degree.

Since I'm new to these concepts, I wouldn't mind if someone points me to some resources or references since it might be hard to provide a succinct answer.

From the article:

  First the image is rotated by an infinitesimal amount $\alpha$. This is done by computing the rotated coordinates of each pixel and interpolating the gray level values at the new coordinates.

Interpolation is key to very small rotations. It allows to sample the image between the pixel center points as if it was a smooth 2-dimensional function. You can choose the interpolation method in skimage.transform.rotate using the parameter order. It controls the order of the interpolation method which (I think) is an interpolating cardinal B-spline. Higher orders give better results but are slower to compute.

If the tangent vector calculated by rotation, subtraction and division is increasingly non-zero as you go further away from the center, then that is an expected result and you can conclude that some rotation took place.


Spectrogram Analysis I have four .wav files and each contain bird chirps.  Given an input .wav file, I need to find out which bird it is.

I have some basic understanding of signal processing. I have used spectrogram to analyze my audio files and now I don't know how to make the frequency spectrum comparison to find out the solution.

Can someone please let me know the process with which I could find out the frequency of sound, bird chirps in my case, programmatically from the spectrogram? I am just a beginner, so I apologize if I said something which wasn't clear

Please check this answer, which describes a few approaches to the same problem. Given that bird song is a monophonic signal (only one fundamental frequency at any point in time - as opposed to polyphonic) - and given that the timbre is irrelevant, the most interesting feature to extract for this classification task is a pitch contour.


I don't understand how the zeros are found of this 2nd order feed-forward filter I realized someone asked a question about the exact same section of a book, which can be found here: Question about zeroes of simple 2nd order FIR filter.
While the answer to this question looks useful to me I have a question about the what the author did to find the zeros. Particularly this part:
$0=1+\alpha_1z^{-1}+\alpha_2z^{-2}$
We know that the zeroes need to be complex conjugates, so can change it to this:
$0=(1-Z_1z^{-1})(1-Z_2z^{-1})$
where $Z_1=Re^{j\theta}$ and $Z_2=Re^{-j\theta}$
The author says that $Z_1$ and $Z_2$ are the zeros, but I don't understand how this can be. Wouldn't $Z_1$ or $Z_2$ have to $=z$ for this to be true?

Below is a polynomial in $z$ :
$$ P(z) = 1+\alpha_1z^{-1}+\alpha_2z^{-2} $$
Zeros of this polynomial is also known as its roots; i.e., those $z_r$ for which $P(z_r)=0$.
Solving for the roots of a general $n$-th order polynomial is not possible. They are numerically found.
However, one alternative way to define a polynomial is to write it as the product of factors of the form $1 - d_k z^{-1}$, where each $d_k$ is a zero/root of $P(z)$; i.e., $(1 - d_k d_k^{-1}) = 1-1 = 0$.
Then your second order polynomial $P(z)$ can be written as
$P(z) = (1-Z_1z^{-1})(1-Z_2z^{-1})$
where the complex valued $Z_1=Re^{j\theta}$ and $Z_2=Re^{-j\theta}$ are the roots of the polynomial. You should solve the second order equation $P(z)=0$ to find out what those zeros are.


Understanding the Threshold Process in Harris Corner and Edge Detector I have read the paper Chris Harris , Mike Stephens - A Combined Corner and Edge Detector about the Harris Corner and Edge Detector, but I didn't understand the point that after low and high threshold they have got 5 regions, specially the point where he has got 2 corner classes, can you please explain that to me with image so I can get that point?

here is the part of the paper which I didn't understand:

  By applying low and high thresholds, edge hysteresis can be carried out, and this can enhance the continuity of edges.

  These classifications thus result in a 5-level image comprising:

    Background.
    Two corner classes.
    Two edge classes.

**Thanks in advance**

After doing the low and high thresholds to the eigenvalues, there will be 2 types of corners:

weak corners which are close to the low thresholding value.
strong corners which are above the high thresholding value.


Implemenation of Length N Moving average filter using block covolution algorithm(overlap and save algorithm) I want to implement a Length N moving average filter using block convolution algorithm(overlap save algorithm)
I am copying block convolution algorithm(overlap save algorithm) from Proakis, DSP using MATLAB 3rd Ed
function [y] = ovrlpsav(x,h,N)
% Overlap-Save method of block convolution
% ----------------------------------------
% [y] = ovrlpsav(x,h,N)
% y = output sequence
% x = input sequence
% h = impulse response
% N = block length
%
Lenx = length(x); M = length(h);
M1 = M-1; L = N-M1;
h = [h zeros(1,N-M)];
%
x = [zeros(1,M1), x, zeros(1,N-1)]; % preappend (M-1) zeros
K = floor((Lenx+M1-1)/(L));         % # of blocks
Y = zeros(K+1,N);
% convolution with succesive blocks
for k=0:K
    xk = x(k*L+1:k*L+N);
    Y(k+1,:) = cconv(xk,h,N);
end

Y = Y(:,M:N)';                      % discard the first (M-1) samples
y = (Y(:))';                        % assemble output

Now using above function i am trying to implement a Length N moving average filter.
My script(m-file) is below
clc
clear all
close all
n=1:5
signal=n*pi;
subplot(221)
stem(n,signal);
xlabel('n');ylabel('amplitude');title('original signal')
noise=rand(1,5);
subplot(222)
stem(n,noise);
xlabel('n');ylabel('amplitude');title('noise signal');title('noise ')
p=ovrlpsav(signal,noise,6),
subplot(223)
stem(p);xlabel('n');ylabel('amplitude');title('corrupted signal')
m=5;
xu=0;
for i=1:m

    x=signal+noise;
    xu=xu+x;
end
xu=xu/m;
subplot(224)
stem(n,xu);xlabel('n');ylabel('amplitude');title('filtered signal')

I am getting a filtered signal similar to input signal as shown by above code,but am i following proper approach?

in psudo pyhon code
const tail = 5

'inside loop where you feed it values

list.push(entry)
list.rotate.left(1)
if len(list)>tail list.pop()
avg = sum(lst) / len(lst)
print avg


Loosing information in decimation - without aliasing We sample each other sample (01010101). Now say we have a signal - that is inside the original sample which has the pattern of (010-1).
Decimation gives us 000 for this sub-signal.
How to prove this is an example of aliasing? Or why else do we have lost information?

The whole point of aliasing and the sampling theoreme is that you cannot (generally) know what an aliased signal was, as you cannot represent infinite bandwidth using finite infornation.

If you have extra knowledge of the input signal (eg that DC is not possible, then you might deduce what a string of 1-1-1 was originally.


How to increase resolution of FFT? I have an input signal in the ( 0 - 20 kHz ) frequency range . When i sample this signal maximum sampling frequency is around 40 kHz . When i calculate the FFT using  1024 points i got resolution of 39 Hz ( resolution = fs/N) . But i need a resolution of 5 Hz . How can i get better resolution?

People interpret the term "resolution" more than one way.  One way is the ability to separate closely spaced frequency lines, and another is to estimate the the "true" frequency of a single tone.  Most DFT interpolators fall into the true value of a tone category.  The problems are not independent of each other but you should understand the problem you are trying to solve.

There is more than one way to interpolate intermediate frequencies from the DFT.

From Peter K's web page, the Fourier Coefficient section links to a survey prepared by Eric Jacobsen.

I've used

Macleod, Malcolm D. "Fast nearly ML estimation of the parameters of real or complex single tones or resolved multiple tones." IEEE Transactions on Signal processing 46.1 (1998): 141-148.

which is a complex interpolator with good result

The basic idea is that there is a largest bin that contains the frequency of interest and the bin leakage is asymmetric. One uses a "fit" to some formula as an interpolation.  THis implies that you need a nominal "high" SNR, and a few bins of separation. Macleod's technique assumes that you used a boxcar window which can be a problem where the noise is far from flat.

You might want to post plots for your problem for further advice.

The other categories  on Peter K's page are also worth looking at.


FFT Matlab - Meaning of Frequency Vector I'm following a tutorial about the FFT. It's well explained but I don't understand the meaning of the frequency vector:

Fs = 150; % Sampling frequency
t = 0:1/Fs:1; % Time vector of 1 second
f = 5; % Create a sine wave of f Hz.
x = sin(2*pi*t*f);
nfft = 1024; % Length of FFT

% Take fft, padding with zeros so that length(X)
is equal to nfft
X = fft(x,nfft);
% FFT is symmetric, throw away second half
X = X(1:nfft/2);
% Take the magnitude of fft of x
mx = abs(X);
% Frequency vector
f = (0:nfft/2-1)*Fs/nfft;

% Generate the plot, title and labels.
figure(1);
plot(t,x);
title('Sine Wave Signal');
xlabel('Time (s)');
ylabel('Amplitude');
figure(2);
plot(f,mx);
title('Power Spectrum of a Sine Wave');
xlabel('Frequency (Hz)');
ylabel('Power');

Can someone explain me what is this for:

f= (0:nfft/2-1)*Fs/nfft;

Here is the link

Specifically,
nfft = 1024
Fs = 150
so f = (0:511)x(150/1024)
ie f will go from 0x(150/1024) = 0 to 511x(150/1024) = 75
So, you are plotting the magnitude of the frequency spectrum from 0 to 75 Hz.
The spectrum is periodic & will repeat for 75 to 150 Hz hence you are plotting from 0 to Fs/2


Undoing a filter operation in MATLAB Out of curiosity, is it possible to actually reverse a filtering on MATLAB. For example, if I have

[b,a] = butter(6,0.6);
dataIn = randn(1000,1);
dataOut = filter(b,a,dataIn);

Is there some way to go from dataOut and obtain dataIn in this example? I was able to more or less get there with the very simple

dataOutprime = filter(a,b,dataOut);

dataOutprime still doesn't exactly overlay on top of dataIn though.

In general it's only possible to implement causal and stable filters. There are exceptions where marginally stable filters are used, but this doesn't apply here. So if you want to invert a given filter, this is only possible if its zeros are inside the unit circle of the complex plane (such filters are called minimum-phase filters). The zeros of the given filter become the poles of the inverse filter, so if the zeros are inside the unit circle, the poles of the inverse filter will also be inside the unit circle, and, consequently, the resulting filter is causal and stable.

The problem with the filter in your example is that it has a sixth order zero at $z=-1$, i.e. on the unit circle. Consequently, this filter is not strictly minimum-phase and it can't be inverted by a causal and stable filter. However, it can be inverted approximately by a causal and stable filter. One way to achieve such an approximate inversion is by minimizing the mean squared error between the ideal and the actual filter response. In practice, it is often acceptable to add some delay to the output, which will make it easier to approximate the desired response by a causal and stable filter.

Note that transients don't play a role. It's the inverse filter's instability that is the problem here.

As an example, take a filter with coefficients

b = [1,0,-0.5]; a = [1,0.3,0.1];

The zeros of the numerator polynomial (with coefficients b) are all inside the unit circle, so the inverse filter is causal and stable.

x = rand(100,1);
y = filter(b,a,x);
x2 = filter(a,b,y);
max(abs(x2-x))    % in the order of 1e-16, i.e., machine precision


STFT - Window Functions I'm looking for your opinions.

I have obtained a Speech Signal and computed a STFT on this signal. I have split the signal into blocks of (100 x 256) and computed the STFT based on different Windowing functions. But, they all bring back similar results (outputted on the same block) I know this might make sense, but, I assumed that because I was using different Windowing functions, it would have a major impact on the results..

Barlett:

Hamming ^^
Note how the block contains no actual speech, just white noise as it's the first block. Does this look normal, or, is there something wrong here?

This is a rather general answer : I would say that, although each window function has specific characteristics that make it more or less suitable for specific situations, in this context you're not going to see a huge difference in the STFT of a windowed speech signal.

In particular, the two windows that you're comparing are quite similar, so I'd expect that you'd get pretty similar STFT results. You might try two windows that are quite different, like a rectangular and a Hamming or something like that, to see just how much difference those make in your STFT, and then use that difference as a sort of upper bound on the difference that you're likely to see from other windows.


Block Diagram for a difference equation I have this little doubt regarding how to draw a block diagram representation of a difference equation. Let us implement $ y(n) = ay(n-3) + by(n-2)+cy(n-1) + x(n)$ in block diagram where $a, b, c$ are constants. What I made is this.

Is this block diagram wrong? They say that you should use only $z^{-1}$ delay multiple times to produce the effect of $z^{-3}$ and $z^{-2}$. But I think that this diagram is also correct. Please help.

That's one way to do it but but typically you would have a single line of three delays in the feedback path and apply $a$ after the first delay, $b$ after the second and $c$ after the third.

And yes, you need arrowheads to make this clear.

By convention the coefficient of the feedback path would be called $a_1$, $a_2$, $a_3$ and the feed forward path would be $b_0$, $b_1$, etc.

See for example https://www.researchgate.net/figure/Direct-Form-I-structure-of-IIR-filter_fig1_236688224


Coherent detection in radar and radio signals Coherent demodulation\detection is usually done via matched filtering in RF signal processing. In radar signals, coherent integration can be done if the phases of the returns are known. In terms of coherent detection in radar signals, does it make any sense if we still use matched filtering, similarly to RF signals? Or, coherent detection means coherent integration in radar signals?

I am new to radar signal processing and trying to clarify some concepts here.

Thanks.

In the context of wireless systems, coherent detection has the definition of detection when the wireless channel is known at the receiver and thus matched filtering is possible and infact optimal at low SNR.

The signals from different antennas at the radar receiver are coherently combined using maximum ratio combining. The MRC is nothing but the matched filter.


Is there any connection between frequency response of a system and spectrum of a signal? For instance, an AR process can be generated through a white noise $w[n]$ passing through an all-poles system (IIR),
$$y[n] = \sum_{k = 1}^{p}a_ky[n-k] + w[n]$$
and we can analyze the IIR system, instead of obtaining the spectrum of the signal, to get some properties of the signal and therefore apply it to some applications, e.g., prediction.

My question is whether all the signals, stationary or non-stationary, can be treated as it is generated via some system driven by white noise, and we can analyze the system, which hopefully can completely characterize the signals.

As mentioned in my comment, look at the Wold Decomposition Theorem :

https://en.wikipedia.org/wiki/Wold%27s_theorem

but actually showing how it relates to your question,

http://www.phdeconomics.sssup.it/documents/Lesson11.pdf.

In simpler terms, your AR system can be represented by a MA system (possibly infinite order) and as the second reference shows, an ARMA system.  The process has more than one representation so there isn't a one to one correspondence between the time series and some system that generated it.

So the answer is yes, there is a correspondence and often the system itself, is of interest, but the time series and system are not uniquely equivalent.

another way of thinking about it is that $H(\omega)$ has a phase and magnitude response that gives $|H(\omega)|^2$ but there are other possible phase responses that results in the same magnitude squared.


What's the theoretical bandwidth? A radio signal is sampled at $f_{s}=5512$ Hz, and $n=2^{14}=16384$ points are used in the FFT to display it as a spectrogram.
The frequency resolution is $\frac{f_{s}}{n}=0.33642578125$ Hz/point, displayed in the spectrogram containing $\frac{n}{2}=8192$ pixels (Nyquist).
Assuming there is no DC offset, the first pixel displayed is $f_{min}=0$ Hz, the last corresponds to $f_{max}=\frac{f_{s}}{2}-\frac{f_{s}}{n}=2755.66357421875$ Hz.
So far, I understand.
But, what's the precise bandwidth?
Is this $\frac{f_{s}}{2}=2756$ Hz, or rather $f_{max}-f_{min}=2755.66357421875$ Hz?

This question seems to highlight some of the perils of translating concepts between the continuous domain and a discrete one.  My understanding is that bandwidth is defined in the continuous domain as the difference between the upper and lower band limits.  Since it is the interval length, whether the interval is open or closed doesn't matter.  Open or closed means if the endpoints are considered part of the interval.  Since we are always saying "up to the Nyquist", considering the endpoints to be open seems logical.

The extra twist in this is your stipulation that the DC component is zero.  This would seem to narrow the bandwidth at the lower end, but to what?  In the continuous case, it would be an infinitesimal amount.

Can this properly be translated to the discrete case you describe?  I don't think so.  There may be a common convention for doing so that I am unaware of.  Perhaps your approach, but then you should strike off the gap at both ends, since you have excluded  DC.

Upon some reflection, I am going to disagree with MBaz again, and say that $f_s/2$ is the theoretical accurate answer based on the continuous definition of bandwidth.

I will agree that these are purely theoretical considerations with negligible practical value, but then again, so are a lot of exam questions.

You should probably strike the word "radio", it stopped RB-J from giving a good answer.  Probably a better answer than what I just gave.

$$ H_k(\omega) = \sum_{n=0}^{N-1}W_N^{-kn}e^{-jn\omega}  $$

$$ W_N^{-kn} = \left( e^{-j \frac{2\pi}{N}} \right) ^{-kn} = e^{j kn \frac{2\pi}{N}} $$

$$
\begin{aligned}
H_k(\omega) &= \sum_{n=0}^{N-1}e^{j kn \frac{2\pi}{N}}e^{-jn\omega} \\
 &= \sum_{n=0}^{N-1} e^{-jn ( -k \frac{2\pi}{N} + \omega) } \\
 &= \sum_{n=0}^{N-1} \left( e^{-j ( -k \frac{2\pi}{N} + \omega )} \right)^n \\
 &= \frac{1 - \left( e^{-j ( -k \frac{2\pi}{N} + \omega )} \right)^N }{ 1 - e^{-j ( -k \frac{2\pi}{N} + \omega )} } \\
 &= \frac{1 - e^{-j ( -k 2\pi + \omega N )}  }{ 1 - e^{-j ( -k \frac{2\pi}{N} + \omega )} } \\
 &= \frac{1 - e^{-j \omega N }  }{ 1 - e^{-j ( -k \frac{2\pi}{N} + \omega )} } \\
 &= \frac{e^{-j \frac{\omega N}{2} }}{e^{-j ( \frac{-k\pi}{N} + \frac{\omega}{2} )}}
 \cdot \frac{ e^{j \frac{\omega N}{2} } - e^{-j \frac{\omega N}{2} }  }{ e^{j ( \frac{-k\pi}{N} + \frac{\omega}{2} )} - e^{-j ( \frac{-k\pi}{N} + \frac{\omega}{2} )} } \\
 &= e^{j \left( \frac{-k\pi}{N} - \frac{\omega (N-1)}{2} \right) }
 \cdot \frac{ \sin \left( \frac{\omega N}{2} \right) }{ \sin \left(  \frac{-k\pi}{N} + \frac{\omega}{2}  \right) } \\
\end{aligned}
$$

For similar math, see my blog article DFT Bin Value Formulas for Pure Complex Tones.  The section called "As N Gets Large" covers the relation between "alias" sinc, aka Dirichlet Kernel, and the normalized sinc function.


Phase Offtset Bandlimited Saw Wave I'm able to make a standard band limited saw wave, but at some point in the past, I remember reading that if you phase offset the sinusoids, you get a different looking wave form that sounds the same.

Trying to implement this by offsetting the phase 270 degrees, I'm seeing a different looking wave form but it doesn't quite sound the same.  It has less volume and sounds a little less bright.

Am I remembering correctly / implementing correctly?  Here's the saw and "other" wave below:

Thanks!

Have you normalized the waveforms differently? The lowest point of the phase-shifted saw should be at amplitude -2 ln(2) / pi = -0.44 if the amplitude of the saw is normalized to -1 .. 1. The positive peak of the phase-shifted saw approaches infinity as you increase the bandwidth. So you could also have a clipping problem. Nonlinearities of the physical world affect the phase-shifted saw more strongly. The two waveforms have the same power spectrum and "complex envelope", but human hearing is more complex and may be able to detect some difference.


Is `xcov()`More Accurate Than `xcorr()` to Compute the Cross Correlation of Signals with Non Zero Mean Values? I have to compute the cross-correlation between two transient signals, which could be decomposed into a trend (with non-zero mean value) + noise.
In particular I am interested in the max of the cross-correlation.
It seems that in MATLAB the cross-correlation is more accurate if it is calculated  using the function xcov() rather then xcorr(), probably as a consequence of the fact that the function xcov() removes the mean value of the signal.
An example of two signals between which I would have to measure similarity can be this one:

Do you know if this makes sense and why?

When trying to measure similarity between signals we're basically building a metric.
When doing so we need to ask what we want to be sensitive about.
For instance, if you don't remove the DC Component (The Mean) and use something like an integral to measure similarity (Correlation / Convolution based) then you are sensitive the added DC component.
For instance, take 2 random white noises and add to each a DC component of 1e6.
Measure their similarity by correlation and compare it to 2 Sine Waves with amplitude of 1e-3 and a relative shift phase of 1e-6 [Rad].
In your eyes, which should be more similar? Which is indeed?
I'd say that in most cases in signal processing, similarity should be measured ignoring the DC component.


Band-edge timing error detector, time accumulator im trying to implement a band-edge FLL and timing locked loop from this paper (https://www.researchgate.net/publication/261499267_Band_edge_filters_perform_non_data-aided_carrier_and_timing_synchronization_of_software_defined_radio_QAM_receivers).
The FLL works fine but i didnt understand the time accumulator and the interpolator.
Usually the interpolator would decimate the signal to 1 smpl/sym.
In this application it is just a variable delay filter?
How do i pick the optimal sample afterwards?

Thanks for your hints and answers. That helped me.
Ive got my code to work now. Its not fine tuned but here is my solution:
-As loop filter i use a PI Filter with a relatively low Ki.
The loop filter is just updated at a 1 Sps/Sym rate when the time accumulator gives an overflow.
-The time accumulator is a modulo interpolation controller
-The interpolator is a farrow structure for cubic interpolation
The IQ plot isn't stationary because the phase isn't locked. But i wont add a following pll anyways because it is a pi/4-DQPSK signal. A bit of phase shifting shouldn't bother me


Ship detection from aerial image I have an aerial image, like this

I want to detect the outline of the ship without shadow. How can I do it? I tried threshold in both graylevel, RGB and HSV. It doesn't work perfectly.

I have used Laplacian of a gaussian filtered image with sigma value of 3.
I have thresholded the LoG image with a high-pass filter.
I have used gray closing morphologic operator.
I have filtered the blobs according to their area.

There are several more ways to do it depending on the images, conditions(sea), illumination, ships etc.
For this image, non-color methods are preferable as shadows are not always predictable.


How to prove these two definitions of the minimum phase transfer function are same? There are so many definitions of the minimum phase transfer function, and these are two of them.

The transfer function of the system which has no zeros or poles at right half plane.
The transfer function which has the minimum phase angle range among the systems which has the same magnitude characteristic.

And these two sentences are describing the same thing. So I want to prove that these statements are equivalent. How to prove it?

PS : Suppose that the system is continuous.

To your second definition it should be added that you only consider causal transfer functions, because it is not difficult to find a smaller phase lag with a non-causal system:

  A minimum-phase system is a causal and stable system with a phase lag that is smaller than the phase lag of any other causal and stable system with the same magnitude response.

Note that a real-valued zero in the left half-plane contributes a phase change of $\pi/2$ to the total phase as we move along the frequency axis from $\omega=0$ to $\omega\rightarrow\infty$:

$$0\le\arg\{j\omega+a\}<\frac{\pi}{2},\qquad a>0,\quad\omega\in[0,\infty)\tag{1}$$

A complex conjugate pair of zeros contribute a phase change of $\pi$.

On the other hand, a real-valued zero in the right half-plane contributes a phase change of $-\pi/2$ as $\omega$ moves from zero to infinity:

$$-\frac{\pi}{2}<\arg\{a-j\omega\}\le 0,\qquad a>0,\quad\omega\in[0,\infty)\tag{2}$$

Note that I've chosen the sign of the term $s\pm a$ such that in both cases the phase is zero for $\omega=0$, which is necessary for a fair comparison between the two cases.

Consequently, exchanging a zero in the left half-plane for a zero in the right half-plane (without changing the magnitude of the frequency response) will always result in an additional phase lag of $\pi$ for $\omega\rightarrow\infty$.

As an example, consider two first-order transfer functions:

$$H_1(s)=\frac{s+2}{s+1}\quad\text{and}\quad H_2(s)=\frac{2-s}{s+1}\tag{3}$$

The signs of $H_1(s)$ and $H_2(s)$ were chosen such that their phases are both zero for $s=0$ (and not $\pm\pi$). The figure below shows the phase plots ($\arg\{H_1(j\omega)\}$ in blue, and $\arg\{H_2(j\omega)\}$ in green):

The pole contributes a phase change of $-\pi/2$ as $\omega$ moves from zero to infinity. The left half-plane zero of $H_1(s)$ contributes a phase change of $\pi/2$, resulting in a net phase change of zero, whereas the right half-plane zero of $H_2(s)$ contributes a phase change of $-\pi/2$, resulting in a total phase change of $-\pi$.

Another way to see the same thing is to note that any causal and stable transfer function can be written as the product of the minimum-phase transfer function with the same magnitude and a causal and stable allpass:

$$H(s)=H_m(s)H_a(s)\tag{4}$$

It can be shown that the phase of a causal and stable allpass is always non-positive for $\omega\in[0,\infty)$, and, consequently, the phase lag of the minimum-phase system is always less than or equal to the phase lag of any other causal and stable system with the same magnitude response.


Simplification after sampling an analog signal $80\textrm{ Hz}$ Sampling rate produces the $x[n] = 3\cos(5\pi n/4)$ signal, the original analog signal is $x(t) = 3\ cos(100\ π t)$.
This is a solved problem used as an example. In the solution it's also mentioned that the resulting signal can be simplified for $|f| ≤ 1/2$ as, $x[n] = 3 cos(5 π n / 4) = 3 cos[(8 – 3) π n / 4] = 3 cos(2 π n – 3 π n / 4) = 3 cos(3 π n / 4)$.

I do not understand how $3cos(2 π n – 3 π n / 4)$ becomes $cos(3 π n / 4)$.
I know that $cos(2πn-3/4πn) = cos(2πn)cos(3/4 πn)+sin(2πn)sin(3/4πn)$.
And $cos(2π) = 1$ , $sin(2π) = 0$ if it wasn't for the "n" variable then yes i think it would equal $cos(3 π/ 4)$.
Why and how does this simplification happen?

Just for explanation: What you experience here, is aliasing. Look at this code and its output:

Fs = 80 # the sampling frequency
f = 50 # the frequency of the cos wave
t = np.arange(0, .1, 1./Fs)  # the sample times
t2 = np.arange(0, .1, 1./(10*Fs))  # higher sampling frequency illustration purposes

plt.stem(t, np.cos(2*np.pi*f*t), label='sampled points')  # plot sampled values
plt.plot(t2, np.cos(2*np.pi*f*t2), label='50Hz cos wave')
plt.plot(t2, np.cos(2*np.pi*(Fs-f)*t2), label='30Hz cos wave')
plt.legend()

What I show is a 50Hz cos wave (blue) and the sampling points when sampling with a frequency of 80Hz (blue dots). I also draw a cos wave with frequency 30Hz on top of the diagram (green).

As you can see, at the sampling points the green and blue line intersect. I.e. when your sampled values do not differ, if you sample a 50Hz or 30Hz wave. This is the aliasing effect: You cannot distinguish between waves that have a different frequency, since the sampling rate $F_s$ is too low. In general, given a sampling frequency $F_s$, waves with frequencie $f$ and $F-f$ create the same sampled values.


Basic questions related to Kalman filter Since I'm a newbie to Kalman filter, so, I had some confusions that I needed to clarify.

What is the difference between measurement and state?
Can they be of different dimensions?

All the examples that I have seen online
     take velocity and position into account while acceleration is
     external control variable. Now the state matrix is multiplied with
      a matrix containing $\Delta t$.

What happens to state matrix and the matrix it is multiplied with to get the updated state if my
  readings are acceleration itself? I am specifically referring to this dataset : Accelerometer dataset
Since readings already have
  noise, would we still add error/noise while calculating the
  updated state?

A state has memory. The value of a state of a causal system reflects the values of all the previous states and inputs. A linear combination of all previous states and inputs, is a linear state equation.

A measurement only reflects the value of the current state and current inputs.  A measurement needs to be noisy. A Kalman Filter estimates the current state from the noisy measurement.

The presence of measurement noise is why it is called a filter. If there is no measurement noise, one uses a state observer, not a Kalman Filter. The KF equations don’t work without measurement noise.

If the measurement noise is a linear combination of the current state and input, which includes noise, the measurement is linear. When the state and measurements are linear, the KF is a linear KF.

Essentially a state has memory and a measurement only depends on the current state. We can write the state evolution as a recursive equation because it is a Markov process.

The dimensional issues are related to what older books called the controlability and observably matricies. More recent books introduce Grammians.

There are variations of the KF that include nonlinear states and  nonlinear measurements. The linear KF satisfies a number of optimality criteria if the matrices and initial conditions are accurate.  The other types aren’t optimal but typically good if the approximations are good.

One does not add error unless you are doing a simulation. you need to know the covariance to actually filter.

The KF can be derived a few different ways, (least squares, orthogonality, bayesian) and the solutions can be written out differently so books often don't exactly agree.  I prefer the innovations form of the KF because the innovations are an online diagnostic.  You can test if you have enough states or too many states.

The questions related to the actual state and measurement equations are really independent of the KF.  A book on linear systems using state variables is where you learn to write those, but many KF books have a chapter or 2 on the topic.  Usually, how to convert from continuous time differential equations to a set of difference equations.  Writing the state equations is really a prerequisite.  Most of the people who ask for help tend to fall in the category of not understanding state variables.

If you have problems understanding your state equations, you should include the specific equations.  No one is watching what you see on the internet. This isn't pornhub.

Kalman Filters tend to work even if you don't know those parameters well.  People tend to spend a lot of time tweaking their filters so expect to tweak a lot.  The theory can be intimidating and one can often build adequate filters without a deep understanding. If it works, it works.  If it doesn't work and you think it should, the theory becomes important.


RMS contrast of image interpretation I have two grayscale images, say A and B, and I want to determine, which one is more contrast. For that I'm calculating an RMS contrast as it is described in Wikipedia. The results are about RMS_A = 0.4 and RMS_B = 0.1. Does it mean that image A is more contrast as its gray pixels are nearer to black/white values? How should I interpreter RMS value to determine contrast?

You defined to compute an RMS. I don't know why you choose it, but OK. How can we probe it as a contrast measure?  If you are considering images with positive pixel values, and consider that multiplying them by a scale factor $a>1$, you extend the range, hence the contrast, or if  $a>1$, you shrink the range, hence the contrast again. How does the RMS behave? A quick computation should give you an answer.

Indeed, those are hypotheses. Honestly, I would not consider RMS as a proper measure, as it does change with a value offset, which to be is substantially iso-constrast.


How to obtain impulse response from the differential equation of a system? This year I'm having trouble with my Signals and Systems class. My major subject is Software Engineering and Electric and Electrical Engineering is my Minor. This question was my previous exam question but I coudn't write anything about it as a solution. If somebody could at least show a way to solve this problem it would be great.

Question:
$$\frac{d^2y(t)}{dt^2}+\frac{6dy(t)}{d(t)}+8y(t)=2x(t)$$
a) Find impulse response $h(t)$
b) Find output $y(t)$ if $x(t)=te^{-2t}u(t)$

I know that $h(t) = \frac {y(t)}{x(t)}$, and i m guessing that i can get that using z-transform but I don't know how to and if it is the correct way.

Thanks in advance...

EDIT:
Here is what I have done so far.

From here, to find the output I think I will use convolution.? y(t) = h(t)*x(t)

I used Laplace transform to find the inverse fourier transform of the function H(jw). What was I going to do if Laplace transform would not be suitable to situation?

It looks like your transfer function is correct, but there's a small mistake in your partial fraction expansion:

$$H(s)=\frac{2}{(s+4)(s+2)}=\frac{1}{s+2}-\frac{1}{s+4}\tag{1}$$

The corresponding impulse response is

$$h(t)=(e^{-2t}-e^{-4t})u(t)\tag{2}$$

The response to $x(t)=te^{-2t}u(t)$ is indeed most easily computed by solving the convolution integral:

$$y(t)=u(t)\int_0^tx(\tau)h(t-\tau)d\tau\tag{3}$$

I leave the exercise of solving $(3)$ up to you, but if I'm not mistaken the result should be

$$y(t)=\frac{1}{4}e^{-2t}\left[2t^2-2t+1-e^{-2t}\right]u(t)\tag{4}$$


IIR - How the gain is calculated? I try to extend my understanding of digital filter design, so I play around with some YouTube videos and the MicroModeler DSP with the following IIR design (based on a video):

I use the following parameters:

Zeros: 0.57 ± 0.78i
Poles: 0.51 ± 0.7i

How does the calculator generate the additional gain with value 0.90134? As I understand it so far the gain value is generated by the equation of the zeros, but I can´t calculate the value for it.

Apparently it's normalized so that the gain at Nyqusit is 0 dB. Potentially they just normalize to the maximum being 0 dB. Gain at Nyquist is given by

$$H(\omega = \pi) = \frac{b_0-b_1+b_2}{a_0-a_1+a_2}$$

whereas gain at DC is

$$H(\omega = 0) = \frac{b_0+b_1+b_2}{a_0+a_1+a_2}$$

Both values are real.


Explanation of OASPL and SEL I have data of different SPL for different 1/3 center frequencies.
Is it correct to sum the pressures at a mean square basis of those SPL values (and then convert them back to SPL) to obtain the OASPL (Overall Sound Pressure Level)?
Furthermore, how to compute a SEL (Sound Exposure Level) from different OASPL over a time period? Just acoustically sum those OASPL values and divide the sum by the number of time steps?
Many thanks in advance!

Yes, you can calculate it by summing all of the 3rd-octave band pressures. Keep in mind, that if the units are pressure then simple sum and conversion to the decibel scale is enough. However, if you have 3rd-octave band levels, then you have to calculate the anti-log of those values, sum them and then convert back to SPL.

Regarding the Sound Exposure Level, you can calculate it from the equivalent SPL $L_{eq}$:

$$SEL = L_{eq} + 10\log_{10}T $$

Where $T$ is the duration of the measurement in seconds. There is a simple intuition behind it. The $L_{eq}$ can be understood as the level of a constant signal, which has the same energy as the varying signal you’ve measured, with the same length. The $SEL$, on the other hand, is the same value, but for a signal that lasts 1 second.

I assume that you will calculate the sum of all pressures (in linear scale) for each of the time instants - you’ve called them OASPL (I don't like this name). Then all it takes to calculate the $SEL$ is to average them, convert to decibel scale and add the duration logarithm of the measurements duration.

Here is a complete example, so you can grasp what I am talking about.
To keep it simple, I will use octave bands instead (8 Hz - 16 kHz), which doesn't actually change anything.

Here are some unweighted levels ($L_{Zeq}$), logged every 0.1 second with total duration being 4.8 seconds.

67.2, 61.3, 46.9, 57.1, 57.7, 51.2, 44.4, 37.7, 35.5, 31.4, 23.2, 19.7
66.7, 60.0, 52.5, 58.0, 56.6, 50.5, 44.8, 38.6, 35.3, 32.7, 24.6, 20.0
64.7, 60.6, 65.8, 56.6, 57.0, 53.3, 45.4, 38.8, 35.0, 34.5, 24.6, 19.4
60.8, 60.5, 66.4, 55.9, 55.7, 52.7, 45.6, 38.7, 33.9, 29.8, 21.7, 19.6
57.3, 63.3, 54.8, 57.4, 58.1, 52.1, 43.8, 39.2, 34.5, 33.5, 26.5, 20.2
53.8, 67.2, 50.7, 56.0, 54.7, 50.4, 43.4, 37.7, 35.5, 35.4, 31.5, 23.8
53.3, 65.4, 51.0, 56.2, 55.0, 51.6, 45.4, 38.6, 35.0, 32.3, 28.5, 21.7
53.6, 61.7, 48.2, 56.1, 56.8, 52.5, 45.0, 39.7, 39.6, 35.5, 27.2, 21.1
52.4, 61.9, 49.7, 56.7, 56.1, 51.5, 45.2, 39.0, 35.5, 31.5, 24.8, 20.5
56.6, 62.2, 44.2, 57.8, 56.7, 50.8, 46.5, 58.0, 51.3, 47.9, 43.4, 40.9
57.6, 63.0, 44.7, 58.1, 55.1, 51.0, 58.3, 48.4, 33.6, 27.8, 22.0, 20.5
60.7, 61.6, 52.7, 56.5, 57.9, 56.7, 46.3, 36.9, 33.3, 29.2, 23.7, 19.9
61.2, 62.3, 51.5, 57.1, 57.2, 51.4, 57.9, 56.7, 49.4, 48.8, 51.0, 53.1
60.4, 59.7, 50.6, 56.7, 55.1, 55.9, 51.7, 38.5, 33.5, 27.5, 21.6, 19.4
60.6, 60.4, 52.8, 55.3, 55.9, 53.2, 43.7, 53.8, 52.0, 50.9, 47.6, 49.8
59.7, 67.1, 48.7, 56.0, 54.9, 51.1, 58.6, 55.0, 41.9, 36.7, 39.2, 28.4
60.5, 68.0, 46.3, 56.6, 55.9, 54.3, 55.9, 55.0, 61.3, 44.5, 43.2, 29.2
61.0, 66.4, 43.3, 57.6, 54.4, 56.9, 52.8, 50.3, 51.1, 36.5, 32.8, 20.5
62.3, 61.2, 44.5, 57.9, 52.3, 57.3, 45.3, 40.0, 43.8, 39.9, 34.5, 28.1
64.7, 59.2, 49.2, 59.0, 58.8, 51.0, 47.5, 52.3, 47.5, 44.9, 38.8, 28.6
66.0, 59.6, 53.2, 62.6, 57.2, 50.6, 46.5, 45.2, 40.7, 35.6, 27.2, 20.0
68.8, 60.3, 51.3, 63.7, 51.7, 48.2, 46.5, 38.6, 34.6, 28.7, 25.3, 21.9
70.2, 60.2, 49.5, 60.2, 54.6, 52.2, 47.8, 38.8, 33.0, 29.3, 29.0, 26.1
69.6, 58.6, 48.9, 60.9, 54.1, 51.0, 45.4, 38.1, 34.3, 31.1, 29.3, 24.1
68.8, 58.7, 47.3, 61.5, 53.2, 48.4, 44.8, 37.6, 35.9, 32.3, 29.2, 25.6
67.4, 61.4, 49.4, 61.3, 55.3, 49.8, 47.2, 40.1, 36.7, 32.8, 26.8, 20.9
68.4, 61.8, 61.3, 60.9, 52.6, 47.5, 45.0, 38.5, 36.3, 32.3, 25.7, 19.8
68.2, 61.9, 64.3, 61.5, 51.5, 49.3, 45.7, 37.8, 33.9, 28.3, 22.6, 19.2
68.9, 59.6, 61.1, 60.4, 54.2, 50.9, 44.2, 39.0, 34.1, 27.9, 22.5, 19.7
69.4, 59.5, 57.3, 59.8, 55.0, 51.4, 45.7, 38.8, 33.6, 29.3, 28.1, 20.7
69.3, 60.5, 54.1, 60.3, 51.7, 50.4, 45.2, 37.9, 34.3, 29.1, 25.4, 20.6
68.8, 59.4, 50.0, 60.3, 50.5, 51.7, 44.7, 36.5, 34.0, 29.3, 25.5, 20.0
67.5, 56.2, 49.4, 60.4, 53.2, 51.0, 47.0, 37.4, 33.5, 28.2, 23.0, 19.5
66.3, 53.8, 51.4, 61.4, 52.6, 50.2, 46.1, 37.0, 33.9, 29.8, 23.5, 20.2
65.1, 58.7, 52.3, 61.1, 51.2, 50.6, 45.6, 38.6, 34.5, 28.8, 21.4, 19.6
67.2, 59.6, 50.0, 60.9, 52.5, 48.4, 45.9, 38.9, 34.3, 31.0, 23.0, 19.9
67.1, 60.3, 50.9, 60.2, 53.4, 52.0, 45.7, 39.2, 35.3, 35.2, 27.2, 22.3
66.9, 61.2, 51.2, 60.4, 52.7, 52.1, 46.7, 38.4, 34.7, 31.1, 24.6, 20.8
65.3, 63.6, 52.1, 59.2, 55.6, 50.0, 46.9, 38.5, 34.4, 35.2, 28.3, 21.2
61.9, 65.7, 52.7, 59.4, 56.6, 50.3, 46.8, 39.7, 53.0, 46.8, 32.8, 26.1
59.7, 65.8, 52.6, 61.1, 54.0, 50.1, 44.1, 39.3, 45.0, 41.1, 28.8, 21.7
57.9, 63.7, 48.4, 60.6, 52.9, 51.5, 46.4, 41.3, 37.6, 35.3, 27.2, 20.9
57.4, 62.6, 48.4, 60.3, 55.2, 52.2, 45.7, 37.3, 34.6, 30.4, 25.4, 20.5
58.3, 61.1, 51.0, 59.9, 53.3, 48.7, 45.7, 37.3, 35.9, 31.5, 29.2, 21.6
60.9, 59.3, 51.4, 59.5, 53.3, 50.5, 46.7, 38.7, 36.3, 31.7, 28.6, 21.5
63.1, 59.7, 48.5, 60.1, 54.9, 51.4, 45.5, 37.8, 37.2, 34.4, 30.3, 22.8
62.5, 57.9, 48.4, 59.9, 54.9, 49.4, 45.1, 39.4, 40.3, 35.4, 30.5, 22.4
60.9, 59.6, 48.8, 59.9, 56.4, 50.7, 43.9, 46.2, 41.0, 39.5, 32.6, 25.0

Since this measurement had been performed using an SLM, I know the exact value of $SEL$, which is $75.3 \mathrm{dB}$.

Let's read it and calculate the instantaneous $L_{Zeq}$ for each frame (the thing that you call OASPL).

#!/usr/bin/env python

import numpy as np

if __name__ == "__main__":
    bands = [8.0 , 16.0, 31.5, 63.0, 125.0, 250.0, 500.0, 1000.0, 2000.0, 4000.0, 8000.0, 16000.0]
    oct_L = np.genfromtxt('octave.csv', delimiter=',')

    oct_P = np.power(10, oct_L/10)   # Anti-log to go back to linear domain
    frame_P = np.sum(oct_P, axis=1)  # Instantaneous pressure in each octave band

    Leq_dt = 10*np.log10(frame_P)    # Instantaneous SPL - OASPL

    Leq = 10*np.log10(np.mean(frame_P))  # Equivalent SPL

    T = oct_L.shape[0] * 0.1
    SEL = Leq + 10*np.log10(T)

    print(Leq_dt)
    print(SEL)

That returns the $SEL$ equal to 75.3 dB and per frame equivalent levels:

[ 69.   68.5  69.6  68.9  66.4  68.1  66.8  64.6  64.6  66.3  66.5  66.4
  67.7  65.4  66.1  69.1  70.3  68.7  66.5  67.7  68.8  70.6  71.2  70.6
  70.   69.5  70.5  71.   70.6  70.6  70.5  70.   68.8  68.   67.6  68.9
  68.9  68.9  68.6  68.6  68.2  66.6  66.1  65.3  65.4  66.6  65.9  65.8]


Need to calculate the VARIANCE of GLCM matrix of an image in MATLAB I have calculated by Matlab the four characteristics of texture derived from the GLCM or co-occurrence matrix (entropy, contrast, correlation, energy) but in what way can the VARIANCE be calculated having the GLCM. Is there some implemented code (matlab function) for that?

This link has details code of how to calculate GLCM variance and other statistical features in matlab.
Reference link https://www.imageeprocessing.com/2019/04/texture-measures-from-glcm-matlab-code.html
For meanX and meanY calculation