Question on Chebyshev filter What does the output signal look like when you insert $2000 Hz$ square pulses to  chebyshev filter with parameters: (Matlab syntax). I tried to do it like this but I don't understand where it goes wrong.

%Making the square pulses

Fs = 1000000;
t = [0:1/Fs:50e-4];
f = 2000;

%square pulse with 50% duty cycle

xt = 1*square(2*pi*f*t, 50);

%Making the filter

Wp = [1000 3000];
Ws = [100 3900];
[n,Wn] = cheb1ord(Wp,Ws,3,20,'s');
[b,a] = cheby1(n,3,Wn,'s');

%Filtering with the coefficient

xt_filtered = filter(b,a,xt);

%Plotting the filtered signal

plot(t,xt_filtered);

But I get something like this which doesn't make sense..?

Allowing for differences between Matlab and Octave ( I only have Gnu Octave to play with,) and correcting the errors noted in the comments, I get the following:
Needed fixes:

Values for Wp and Ws divided by half of sampling rate - the cheb1ord function requires frequencies be given relative to the Nyquist frequency.
Removed 's' (not allowed for Octave anyway.)

To make the code work for Octave, I had to make the following changes:

Remove duty cycle value from square
Remove 's' from cheb1ord and cheby1

Code for Octave (should run under Matlab as well:)

%Making the square pulses
Fs = 1000000;
nyq=Fs/2;
t = 0:1/Fs:50e-4;
f = 2000;

%square pulse with 50% duty cycle
xt = square(2*pi*f*t);

%Making the filter
Wp = [1000/nyq 3000/nyq];
Ws = [100/nyq 3900/nyq];
[n,Wn] = cheb1ord(Wp,Ws,3,20);
[b,a] = cheby1(n,3,Wn,'s');

%Filtering with the coefficients
xt_filtered = filter(b,a,xt);

%Plotting the filtered signal
subplot (211);
plot(t,xt);
subplot(212);
plot(t,xt_filtered);


How we define frequency in an image? How we define frequency in an image?

Does the repetition of a single gray level defines the frequency?(Example: can a black(single color/constant region) background image be called high frequency image?) I think it is not so.

Or

Does the frequency depends on number of changes in the gray levels in an image? If so Why? I couldn't able to relate the image frequency with the signal frequency where we say high frequency signal when the same cycle is repeated many times in a unit time.

If an image has large values at high frequency components then the data is changing rapidly on a short distance scale. e.g. a page of text

If the image has large low frequency components then the large scale features of the picture are more important. e.g. a single fairly simple object which occupies most of the image.

Thus low frequency correspond to slow varying information ( continuous surface)
high frequency correspond to quickly varying information ( eg. edges)


Detect lines transitioning to blobs in images I have the following images (in high resolution, uncompressed). When the process starts, the left image is displayed and it ends something similar to the right image.

The first image shows labyrinth-kind-of strucutres with sometimes only bright borders, sometimes dark borders and sometimes both. The second, middle image shows the transition to the last image: The lines are getting shorter and disconnect. In the end nearly all of the lines will be gone and little circles will be present.

My goal is to detect the transition from lines to blobs, so basically the center image. I need to know when the lines are vanishing or when the blobs are appearing (let's say there are more than 5). The detection doesn't have to be super precise, I just want roughly to know when the first bubbles start coming. Also the detection should be as fast as possible.

What I have tried so far

I am very new to signal processing, so I'm not very sure how to start. I am using python.

My own solution (doesn't work)

In my naive way I thought "These contrasts should be easy to detect.". So I played around with blurring the image a little bit (to remove noise, even though there is very few noise) and then comparing each pixel value to a threshold. I encountered three problems:

There is a (varying) gradient in the background.
The contrast of the top left lines compared to their environment is a different than the contrast
of the lines on the bottom right (because of the gradient)
The lines are bright or dark

I then tried to somehow get the "local" contrast, so to compare the pixel value to the pixels around. I created a mask that selects only parts of the image, but this was (in my implementation) super slow. So I stopped with this way.

In addition I think it is harder to detect vanishing lines than appearing bubbles (?).

Blob detection

I found a blob detection page on scikit-image.org. I tried to use their algorithm and it is pretty good. I get the following result:

I used the transition image (the middle one of the top three images). The Difference Of Gaussian seems like what I want. But it is still detecting bigger "blobs" as you can see on the bottom right of the image wich I don't want. Also I'm not sure if this is the "default way to go" or if there are a lot better things which I'm just not thinking of.

Edit: The blob detection is very slow on my original data. So this is not an option except it can be speeded up a lot.

Canny edge detector

I tried to use the canny edge detector method as mentioned in the comments. But after some trying around I didn't get any good results. I end up with either way too much edges or no edges. Even after playing around and trying to smooth the image before or subtract, divide or multiply a smoothed image. But that may also just be because I don't know exactly how to do this. Out of my very limited knowledge I would say that the contrast is not good enough. But I don't know how to fix that.

So, to conclude: What is the best method to detect those blobs? And maybe, in addition: Is there a best python module to use? And if so, which one is it?

If you are not interested in picking out the frame that bubbles start forming and since the transition is (clearly) from order (clearly formed shapes) to disorder (mostly noise), you can use a feature based approach. That is, characterise the whole image with a few numbers.
An obvious feature here is image entropy. The typical (but not the only) way to calculate entropy of an image / signal is through its normalised histogram (you need to normalise it, to get to a probability).
Are there any differences between the histograms of the images?
Yes there are (at least, the ones that come with the question):

From left to right, we have the image with the well defined shapes (ordered), the middle state and on the right the final state where there is mostly background plus some noise.
The striking feature here is that when the image is ordered, the ridges that are well defined because they create contrast with their background create a bimodal histogram. One mode corresponds to the background (the lower one) and one mode corresponds to the "foreground" (the one that forms a "shelf" towards the right) or the worm-like ridges.
You can catch these differences in many different ways including (of course) distance from Gaussian distribution, skewness, etc. Where the image is going to end up is a set of pixels distributed normally, any deviation from that "profile" must be due to formation of those ridges (in this particular context).
If you try to measure (Shannon) entropy over these images, you get the following three values:
    Ordered: 7.2331079063161425
    Middle:  6.734680695379271
    Noise:   6.591988661960253

(These measurements come from skimage.measure.shannon_entropy() )
This number looks at the frequencies of occurence of the pixels and returns the number of bits required to describe the value of a pixel.
You might feel here that these numbers are contradicting: Low entropy means more order, high entropy means less order. "Less order" here means closer to an absolutely uniform distribution. That is, you cannot really predict what comes next, the image is entirely random.
But at the same time, this description might also help understanding what is going on here: The "Ordered" image alternates between two means (more preference on the lower ofcourse) and spans a wider range of highly probable grey tones and the "Noisy" image is basically $128 \pm \approx 70$ with diminishing probability.
Shannon entropy evaluated in this way looks at the pixel values independently of each other. But they are not. The meandering structure forms straight and curved points to the extent that if one pixel is "on", it is highly likely that the pixel next to it is on as well. So, we are not talking about single pixel probablities any more.
To "catch" this kind of differences, you need a measure of entropy that is sensitive to them. One of them is Entropy as calculated through the co-occurence matrix.
The name is kind of self-descriptive but very briefly: The co-occurence matrix describes how the values of two different pixels that are some distance appart co-vary. Additionally, since we now talk about the distance between two pixels, we can also define a bearing. This offset, implies that the co-occurence matrix describes co-occurence along a direction.
The corresponding values now are:
    Ordered: 14.197350713914734
    Middle : 13.128723041575485
    Noisy  : 12.824963198002461

(These measurements come from mahotas.features.haralick() with return_mean=True. Otherwise, you get the value of each feature for each direction considered (typically, the 4 directions around a pixel))
We observe a similar trend.
The difference between the two ways of computing the entropy here is robustness and sensitivity. We do not have a long enough sample here, but it is likely that the haralick features would be more sensitive to structural changes than a single histogram.
Therefore, you can get a long enough set of images from the experiments, "map" how the value of this feature varies and then set up "alarms" by setting thresholds around specific ranges.
Hope this helps.


Difference of Gaussians I use the following formula for DoG:

$$\frac{1}{\sigma} (\frac{x^2}{2\sigma^2}-1.0) e^{\frac{-x^2}{2\sigma^2}}$$

What is the relationship between this formula and the difference of two Gaussian functions?

Can someone show me how to apply DoG on a 1D signal in MATLAB? My input can be considered as a 1D real-valued function, and I applied the above formula to detect blobs. The results make sense except that I don't know why I get many local minima in a relatively flat region (i.e., function values at points in the region almost stay constant). Since I don't have enough signal processing/image processing background, I want to see what a constant signal look like after being convolved with a DoG filter.

The difference of gaussian (DOG) is the convolution of input image by difference of two gaussians usually with different standard devitations($\sigma$). The basic idea behind this is to capture edges or gradients in the images that are simplified by the gaussian with larger $\sigma$  but preserved by the smaller gaussian. This is form of scale space derivative, which basically analyses different frequency components in the image, by changing the $\sigma$s  of gaussians. Thus, its like a band-pass filter that rejects certain frequencies in the input.

If we consider an constant image with its value centered around the means of the gaussians, the result of the DOG would be a zero - indicating no gradient.Try this out for a basic tutorial.


Bias of random variable I have a really simple question of definition.

Consider that I have an estimator $\widehat{a}$ and a variable $a$.

I know that we define the bias of $\widehat{a}$ by $B(\widehat{a})=\mathbb{E}(\widehat{a})-a$.

I just want to know in the case where $a$ is a random variable, if the definition is this one:

$B(\widehat{a})=\mathbb{E}(\widehat{a})-\mathbb{E}(a)$.

In the case $a$ is a random variable (r.v.) we distinguish between two types of unbiased: (hereinafter I'm switching to $\theta$ and $\hat{\theta}$ instead of $a$ and $\hat{a}$ respectively):

Unbiased in the strong sense: given by $E\left \{ \left (  \hat{\theta} - \theta \right )| \theta\right \} = 0, \quad \forall \theta$. Most often tha n not it is hard to satisfy this definition.
Unbiased  in the weak sense: given by $E\left \{  \hat{\theta} - \theta \right \} = 0, $ where the expectation is executed on all random terms inside. The meaning of this criterion is that you want the estimator to be unbiased "in the mean" and not for all values of $\theta$. This criterion can be calculated using the law of total expectation (smoothing theorem).


Derivation of Yule-Walker Equation for LPC I have derived the Yule-Walker equation as shown below.

\begin{align}
\hat s(n) &= -\sum_{k=1} ^ p a_k s(n-k)\\
e(n) &= s(n) - \hat s(n) = s(n) + \sum_{k=1}^pa_k s(n-k)
\end{align}

In order to minimize the power of the error: $\displaystyle
\frac{\partial E\left[\lvert e(n)\rvert^2\right]}{\partial{a_k}} = 0 \quad \text{and}\quad \lvert e(n)\rvert^2 = e(n)e^*(n)
$

$$
E\left[\left(s(n) + \sum_{k=1}^p a_k s(n-k)\right)e^*(n)\right]= E\left[s(n)e^*(n)\right] + \sum_{k=1}^p a_k E\left[s(n-k)e^*(n)\right]
$$

\begin{align}
&\frac{\partial E\left[\lvert e(n)\rvert^2\right]}{\partial{a_k}} = E\left[s(n-k)e^*(n)\right]=0\\
&E\left[s(n-k)\left(s^*(n)+\sum_{i=1}^p a^*_is^*(n-i)\right)\right]=0\\
&E\left[s(n-k)s^*(n)\right]+ \sum_{i=1}^p a^*_i E\left[s(n-k)s^*(n-i)\right]=0\\
&\sum_{i=1}^pa^*_i r(i-k) = -r(k)
\end{align}

A few things confuse my mind. So I have 4 questions:

Is this derivation true?
I am not sure why we do the partial derivation in terms of $a_k$.
$e^*(n)$ includes $a^*_i$ coefficients, too. Why is $e^*(n)$ not affected by the derivation operation?
What is the meaning of all these complex conjugates here? How do they affect the numeric calculation?

Yes, your derivation is correct.
There is just a minor inconsistency in your definition of the auto-correlation function $r(k)$. If you choose to define
$$r(k) = E[s(n-k)s^*(n)]$$
then $E[s(n-k)s^*(n-i)]$ equals $r(k-i)$ or, equivalently, $r^*(i-k)$, but not $r(i-k)$. So the last equation should read
$$\sum_{i=1}^pa^*_i r^*(i-k) = -r(k)$$
You take the derivative w.r.t. to $a_k$ because these are the coefficients that you can choose as to minimize the error. They are the unknowns in your problem.
The fact that $e^*(n)$ can be viewed as a constant when taking the derivative w.r.t. $a_k$ has to do with how computing the gradient w.r.t. complex variables works. It is equivalent to taking the derivatives w.r.t. to the real and imaginary part of the coefficients separately, but the complex technique is more elegant and easier (at least as soon as you understand it). I've explained this technique in this answer.
The complex conjugates just pop up because you consider $|e(n)|^2$ in your optimization objective: $|e(n)|^2=e(n)e^*(n)$. You wouldn't need them at all if you only considered real-valued coefficients $a_k$ and real-valued signals $s(n)$.


Convert Triangle Wave to Sine Wave I have this triangle wave as input:

and I want to derive a sine wave out of it. How do I go about it?

PS: I am doing this in Lua, so any piece of code close to that language would be appreciated.

Assuming $x(t)$ swings between -1 and 1, $\sin \frac{\pi}{2} x(t)$ might do the trick.


OpenCV warpPerspective implementation I've implemented the least squares method to find the homomorphic image to fix the rotation and projection in an image.

Now I'm trying to implement the OpenCV warpPerspective method in order to "fix" my image, my python implementation is like:

def fix_image(img, t):
    new_image_map = {}
    minx, miny = img.shape[0], img.shape[1]
    maxx, maxy = 0, 0
    for i in range(img.shape[0]):
        for j in range(img.shape[1]):
            xy = np.array([i, j, 1], np.float64)
            uv = np.matmul(t, xy)
            uv = uv / uv[2]
            minx = min(minx, uv[0])
            maxx = max(maxx, uv[0])
            miny = min(miny, uv[1])
            maxy = max(maxy, uv[1])

            new_image_map[int(uv[0]), int(uv[1])] = (i, j)

    minx, miny = int(minx), int(miny)
    maxx, maxy = int(maxx), int(maxy)
    final_img = np.zeros((maxx - minx + 1, maxy - miny + 1)) \
        if len(img.shape) == 2 else np.zeros((maxx - minx + 1, maxy - miny + 1, img.shape[2]))

    for k, v in new_image_map.items():
        final_img[k[0] - minx, k[1] - miny] = img[v]

    return final_img

I know that I still need to interpolate the empty points but the problem is that the shape is not right, I'm checking the results by comparing with the actual OpenCV implementation.

dst = cv2.warpPerspective(storm_img, tr, (1448, 1456))

As you can see, it is far from the expected.

resp = fix_image(storm_img, transformation)

I do not want to just use the OpenCV method because I want to learn how to implement it.
What am I getting wrong?

The problem in your implementation is that it returns the tensor of float, while an image must be a tensor of int. Because of that, your rendering library, which I assume is matplotlib, cannot correctly plot an image.

To fix that you need to specify the type of final_img explicitly. That is, you need to add a parameter dtype=np.int when creating final_img, giving us the following:

final_img = np.zeros((maxx - minx + 1, maxy - miny + 1), dtype=np.int) \
        if len(img.shape) == 2 else np.zeros((maxx - minx + 1, maxy - miny + 1, img.shape[2]), dtype=np.int)

You can also refer to this python implementation of cv2.warpPerspective for more details.


How to calculate CWT shannon entropy? I am trying to calculate shannon entropy of CWT. I am not sure if I am doing it right.
Assume that $W(a_i,t), i=1;2;...;M$ is a set of wavelet coefficients. The Shannon wavelet entropy is calculated by:
$E=-\sum_{i=1}^{M}d_i log(d_i)$ $\rightarrow$   where $d_i=\frac{|W(a_i,t)|}{\sum_{j=1}^{M}W(a_j,t)}$
I am confused how to calculate $E$. for example I have a coefficient matrix with size of $M\times N$, $M$ is scales number and $N$  is time segments. first I have to calulate $d_i$, this is my main problem. this is the wavelet coefficient matrix :
 $W_{M\times N} =
 \begin{pmatrix}
  w_{a_1,1} & w_{a_1,2} & \cdots & w_{a_1,N} \\
  w_{a_2,1} & w_{a_2,2} & \cdots & w_{a_2,N} \\
  \vdots  & \vdots  & \ddots & \vdots  \\
  w_{a_M,1} & w_{a_M,2} & \cdots & w_{a_M,N}
 \end{pmatrix}$

hmm i am pretty sure i am wrong, can anyone help?
for example tell me how can i calculate $d_4$?
Here I have right a little Matlab script to calculate shannon entropy of CWT.
Is it right or wrong? and what should I do?

 [M,N]=size(coeffs);
for js=1:M
Ej(js)=sum(abs(coeffs(js,:)));
end;
Etot=sum(Ej);
Pj=Ej./Etot;
%shannon entropy
shan_entr=-sum(Pj.*log(Pj));

Yes. Your code does compute the same thing as the formula from the paper.


How to "rotate the phase"? I have problem, not sure how to call it but I think it's something called "phase rotating", but not sure.

I make fast Fourier transform on simply square wave, and then I make inverse fast Fourier transform - just for fun :)
and after this forward/backward FFT on square wave I expect on a wave graph something like that:

But I get that:

Could anyone help me, which parameter I should change, to fix it?

Ok I give the answer by myself, but it’s only half of answer. As I know we can rotate $ 90^o $ the phase by multiplying the wave by $ j $ - imaginary number. And to rotate it to $ 180^o $ we need multiply by $ j^2 $ (or just $ -1 $ which is the same). To rotate $ 270^o $ we need to multiply by $ j^3 $, and so on...

But if we want rotate something between? For example $ 45 ^o $ or $ 30^o $. What we need to multiply by?


Sample and Hold Circuit at Digital to Analog Converters Does DACs contain sample-and-hold circuits? I am a bit confused
thanks

The output of a DAC represents a zero-order hold (equivalently called a sample and hold) of the digital input samples: each sample is held at the value in the analog domain for the time duration between samples. There are also other DAC implementations such as a RTZ (Return-to-zero) DAC which brings the signal back to zero midway in the duration.

The zero-order hold is equivalent to convolving in the time domain the digital samples, which represent impulses in the time domain, with a rectangular pulse as wide as the time duration between samples. Convolution in the time domain is multiplication in the frequency domain, and therefore the DAC output will have a "droop" over frequency with the magnitude going to zero as the frequency approaches the sampling rate. The RTZ implementation moves the first zero out to twice the sampling rate, at the expense of signal loss.

See pages 13 and 15 of this TI presentation for more details: https://training.ti.com/sites/default/files/docs/TIPL4705%20-%20DAC%20Output%20Response.pdf


What happens on signal during aliasing? So I have signal with frequency of $15Hz$ and sample rate frequency is $20Hz$.

What happens if we sample the signal at a frequency that is lower that the Nyquist rate?  We will have aliasing. Sample rate frequency should be at least bigger than $2*f_{max}$.

I understand that part but I'm confused how that will look like. Can someone help me with that?

You have aliasing for sampling signals of any frequency: a signal of 5Hz sampled at 20Hz is indistinguishable from a signal of 15Hz of opposite phase.

It's common in audio usage to work with "lowpass signals" only and consider "aliasing" for frequencies more distant from DC, but in reality sampling of anything is ambiguous by multiples of the sample frequency regardless of the original signal.  You can perfectly well sample signals of 15Hz with a sampling frequency of 20Hz as long as you make sure that the signal only contains frequencies between 10Hz and 20Hz (for example).

For any detected signal component, increasing/decreasing its frequency by multiples of the sampling frequency will not cause any change at the sampled time points.  Hence you want to prefilter the signal such that adding multiples of the sampling frequency to any signal component leads to something that cannot be present in the signal: that way the sampled values represent the signal uniquely.


How to decide filter order in Linear Prediction Coefficients ( LPC ) while calculating formant frequency features? I am new in signal processing and trying to calculate formant frequency features for different .wav files.

For calcuating formant frequency, I need three parameters values :

Linear Prediction Coefficients ( LPC )
root
angle

I am trying to calculate Linear Prediction Coefficients ( LPC ) using librosa.core.lpc in python. it takes two parameters:

librosa.core.lpc(y, order)

I have Y but I don't know how to calculate order, I have many .wav files and I have to set order to extract features from all the files. How to determine the order for all wav files to calculate LPC?

Next two things root and angle can be calculated easily like this :

 rts = numpy.roots(A)
 rts = [r for r in rts if numpy.imag(r) >= 0]

 angz = numpy.arctan2(numpy.imag(rts), numpy.real(rts))

 # Get frequencies.
 Fs = spf.getframerate()
 frqs = sorted(angz * (Fs / (2 * math.pi)))

But I am already stuck in the first step.

Thank you in advance!

You have to define the order based on your application, usually for speech you'll use and order of 10, this definws the order of your all pole filter and it's related with the amount of formants tha you will see in the impulse response of the filter.

https://www.dsprelated.com/showthread/speechcoding/965-1.php


problems calculating phase propagation (in C) I'm coding a phase vocoder in C for time-stretching audio signals. The only piece that is currently not working is the phase propagation.

I take overlapping frames from an input signal, apply a window function, and perform FFT.
Then i calculate the phase propagation (the one I'm stuck with)
Finally I perform iFFT and overlap-add the output frames with a different hop size.

If I skip step 2, the phase propagation, the output sound is correctly time-stretched but I can hear a phaser effect. My assumption is that the missing piece will correct this. However when I apply my current implementation the sound is terribly distorted.

I appreciate if anyone can give some good pointers on how the phase correction can be implemented.

It might be too much to ask, but this is my current implementation should you see something weird.

// _prevPhaseX[] = FFT phases from previous frame (n-1) (propagation input)
// _prevPhaseY[] = FFT phases from previous frame (n-1) (propagation output)
// _phase[] has the FFT phase information from current frame
static int first = 1;
if (first)
{
    first = 0;
    for (int k=0; k<half; k++)          // half = window length/2 = number of bins
    {
        _prevPhaseX[k] = _phase[k];
        _prevPhaseY[k] = _phase[k];
    }
}
else
{
    // int length = window length - e.g. 1024
    // float Ra = analysis hop size - e.g. 256
    // float Rs = synthesis hop size - e.g. 512
    // float PI = PI
    // float DPI = 2*PI
    for (int k=1; k<half; k++)
    {
        // compute Delta Phi(k)
        float phase = _phase[k];
        float omega = DPI*(float)k/(float)length;  // length vs length/2?
        float phasediff = phase - _prevPhaseX[k] - Ra*omega;
        _prevPhaseX[k] = phase;   // Save value for next frame

        // get principal argument
        while (phasediff>PI) phasediff -=DPI;
        while (phasediff<-PI) phasediff +=DPI;

        float freq = omega + (phasediff)/Ra;
        _phase[k] = _prevPhaseY[k] + Rs*freq;
        _prevPhaseY[k] = _phase[k];   // Save value for next frame
    }
}
// the _phase[] vector (and magnitude[]) goes to IFFT (after conv. to rectangluar)

Thanks!

you don't wanna do this on a bin-by-bin basis.  here's a paper by Miller Puckette that is seminal regarding this issue.  what you want to do is identify specific frequency components (perhaps each lobe in the spectrum would be a specific frequency component) and multiply that entire $n$th lobe by $e^{j \phi_n}$  where $\phi_n = \omega_n \ T_n$  and  $\omega_n$ is the frequency of the $n$th component and $T_n$ is the time adjustment to that component to match its phase to the corresponding component in the previous frame.

each bin gets multiplied by $e^{j \phi_n}$, but for most bins that factor is the same as for an adjacent bin.


Phase of signals in real and complex DFT I have been reading chapter 8 - The Discrete Fourier Transform of the book The Scientist and Engineer's Guide to Digital Signal Processing and chapter 12 - Discrete Fourier transform of the book Continuous and Discrete Time Signals and Systems. I'm not sure if I have understood real DFT properly.
From what I understood, in real DFT, if the samples in time domain have size $N$, then the output is consist of $\frac{N}{2} + 1$ complex number. Denote the output as X[ ], then Re X[ ] (real part of the complex numbers) is a list of amplitudes of cosine wave, and Im X[ ] (imaginary part of the complex numbers) is a list of amplitudes of sine wave. As is illustrated in the image:

And all of these waves have phase equal to $0$.
Yet if we use the forward DFT formula $$ X[r] = {\sum}_{k=0}^{N-1} x[k]e^{-j(2 \pi k r / M)} $$ to compute ($M$ is usually considered to be equal to $N$), the output is a sequence consisting of $N$ complex numbers, for a input sample with size $N$. And each of those complex numbers represents a frequency component. Denote a complex number in the list as $a + bi$.
Is it the case that the frequency components have amplitude $A = \sqrt{a^2 + b^2}$, phase $\varphi = \arccos{(\frac{a}{A})}$ and they are shifted sine waves?

And all of these waves have phase equal to 0.

No. The phase will be zero for $X[0]$ (DC -zero frequency component) but the phase is given by $arctan(Im[i]/Rx[i])$, which usually doesn't result in zero.

Yet if we use the forward DFT formula
$$X[r] = {\sum}_{k=0}^{N-1} x[k]e^{-j(2 \pi k r / M)} $$
to compute ($M$ is usually considered to be equal to $N$), the output
is a sequence consisting of $N$ complex numbers, for a input sample
with size $N$.

N output complex numbers requires just as many input complex numbers (which means you need N real numbers and N imaginary numbers (all zeros) as input, which together make up the complex input).
If, the input is ,however, purely real (N real numbers), and you don't provide an imaginary part, your real fft algorithm will simply treat your input as complex (for instance, the even numbers will be treated as real and odd as imaginary) but now you'll get only $N/2+1$ complex numbers at the output.
I

s it the case that the frequency components have amplitude $A =
> \sqrt{a^2 + b^2}$, phase $\varphi = \arccos{(\frac{a}{A})}$ and they
are shifted sine waves?

Yes, but the the phase fomula is actually $\varphi = \arctan{(\frac{b}{a})}$
EDIT: As mentioned in a comment below, one generally needs to use $atan2$.


Estimating the SNR of 2-FSK I have some iq data recordings of a 2-FSK waveform that I've created periodograms from. How do you estimate the SNR of this waveform in the frequency domain? I've heard two opinions:
One:

Integrate the signal band power over the excursion
Find a clear area of the spectrum and integrate the noise power over the same bandwidth
Divide the two quantities: SNR = signal power/noise power
This seems weird because FSK is just two tones and integrating the signal over the excursion is just integrating a lot of noise and very little signal

Two:

Find the peak signal power of one of the tones
Estimate the median of the noise floor power
Divide peak signal by median of noise floor: SNR = peak / median
This also seems weird because depending on the snapshot of data, the amplitude of the tones is different. SNR isn't changing but this method doesn't give a constant SNR.

If I could collect the same iq data with the waveform turned off (just receiver noise), would it make the calculation any simpler?

I've used the first technique before. The idea is that you have some chunk of time domain samples where you know only noise is present and use that to get an estimate of the noise power, $P_n$. Then you determine the total received power, which will have contributions from both signal and noise, $P_{s+n}$. When you take the ratio of these, you get:
$$\frac{P_{s+n}}{P_n}=\frac{P_s}{P_n}+1$$
For high SNR, the extra $1$ doesn't matter too much but to be exact the SNR estimate should be:
$$\hat{\text{SNR}}=\frac{P_{s+n}}{P_n}-1$$


How can we combine 2D image with 3D image for the same sense? We have 3D point cloud data captured by a laser scanner and need to register it with 2D images to get further information. In our case, the 2D image and 3D points are not overlapping and therefore we need to preprocess them so as to extrinsically calibrate the 3D scanner with the 2D image sensor (camera).

For the same modality, such as 2D-2D, this is done with good results using standard techniques, but for different data space, we still face challenges in this task.

These kind of images can be used in auto-driving cars and other applications.

Formally speaking, you would like to extrinsically calibrate the laser scanner to the 2D image. I have taken the liberty to edit your question to reflect that.

Here is how my initial approach would be:

Calibrate the intrinsics of the 2D camera. For that, just use OpenCV. You should store the intrinsic parameters: focal lengths, principal point and distortion coefficients.
Take an object which you could easily locate in 3D and 2D. Such an object could be any 3D model for instance if you can perform pose estimation both in 3D and 2D. Because it would be difficult to do this for free-form objects, I would propose to use a sphere-object which appears to be ellipse-like in the camera, and a sphere in the laser scanner. For instance, take an unused pilates-ball.
Place the pilates ball in the overlapping area of the point cloud and the 2D image. Then capture multitudes of images (at least about 30) by moving the ball in the overlapping region to various locations. At each location, make sure that the ball occupies about $1/4^{th}$ of the image and record the 2D and 3D frames.
Detect the ellipses in all the images, i.e. for each image, one ellipse. Use any open source tool, for instance 1D-hough transform to do that. I strongly suggest to refine the ellipse to subpixel accuracy, for instance via this method (source code). Record the center of the sphere (ellipse) in 2D $\{(x_i,y_i)\}$.
Detect the spheres in all of the laser point clouds. Again, many packages do exist. Even MATLAB has one method. Record the detected centers of the 3D sphere $\{(X_i,Y_i,Z_i)\}$. Again, a non-linear refinement of the sphere in 3D can help you in getting more accuracy. Check this one out for a complete pipeline on this if you like more sophisticated approaches.
Now you have an easy problem to solve: You have established 2D-3D coorespondences in the form $\{(x_i,y_i)\}\sim\{(X_i,Y_i,Z_i)\}$, have the calibration information of 2D camera and thus, all you need to do is to solve the PnP problem. Again many open and closed source software exists such as OpenCV's solvePnP or MATLAB's estimateWorldCameraPose.

The result will be a rotation matrix $R$ and a translation vector $t$ that maps the 3D scanner's coordinate frame to the image domain, thereby solving the registration problem.

I hope instructions are clear and this helps.

Note that, many things can be done to make this procedure easier, more robust, more accurate and etc. This created quite a big interest from academia, and if you dig a bit, you will come across a lot of sources. Nevertheless, for the beginning, the aforementioned procedure is pretty easy to implement.


how to guess what the interference source it could be At work, we are in an industrial area, and there is always a signal that I receive inside my office which interferences with my tests.

This signal is in the 433 MHz and when I see the spectrum, there are 2 peaks around this center frequency which are shown sort of constantly pulsing, like if someone were pressing a button every 5 seconds (but this is a 24h interference, so it must be something else).

Also, because of the strength of this signal, I wonder if it's something inside my office but we only have routers, wireless phones...I can't think of any wireless device inside my office that could be transmitting at this frequency.

I only can thing about the garage door openers but, as this interference is constantly there, it must be something else.

What sort of device could be doing this?

I read about it and I found out that it could be our alarm system because the power is high and the transmission is constant so I run the technical support and they confirm, so I guess, problem solved.

However, if someone else had the same issue and found that it was something else, please share :) I read it could be also the thermostat but couldn't find the specs of ours and unless it transmits also when it's not working, I don't think it could be that!


Confusion Related to the Laplacian Filter I have this confusion related to Laplacian filter. It states that it uses second derivative and finds the point of zero crossing(where the second derivative is zero) for which the first derivative has the high peak. I have seen them using kernels like

0 1 0
1 -4 1
0 1  0

I am not sure how this kernel finds the point of zero crossing. I mean definitely this kernel is equivalent to the second derivative.

Suppose I have points something like this

0 255 255
0 255 255
0 255 255

Definitely, it has an edge at the middle. But if I convolve the kernel with this image for example. I won't get zero result. So how come its called zero crossing

The Laplace filter is the sum of the differences of neighbourhood pixels compared to the central pixel. Its output will be zero in regions where the intensity is constant, and non-zero when there is a transition in intensity (i.e. an edge).

For simplicity's sake, consider what happens in one axis only, where the filter kernel would be [1 -2 1]. If we have a sequence of values with a step transition then the input and output will look like this:

input:  10   10   10   10   10  100  100  100  100  100

output:       0    0    0   90  -90    0    0    0


Low Pass Filter to Maintain Edge Information I am looking for a kernel acting as a low pass filter that satisfies these conditions:

$$K(-\mathbf{u})=K(\mathbf{u}) \tag{1}$$
$$K(\mathbf{u}) \ge K(\mathbf{v}), \;\; \text{if} \;\; |\mathbf{u}|<|\mathbf{v}|, \;\; \text{and} \;\; \lim_{|\mathbf{u}|\rightarrow\infty} K(\mathbf{u})=0  \tag{2}$$
$$\int K(\mathbf{x})d\mathbf{x} = 1 \tag{3}$$

In the reference paper, the author suggested a Gaussian kernel that is:

$$K_\sigma(\mathbf{u})=\frac{1}{\left(\sqrt{2\pi}\sigma\right)^n}e^{-|\mathbf{u}|^2/2\sigma^2} $$

with a scale parameter $\sigma > 0$.

The Gaussian kernel is very good at approximating the required properties. But the kernel reduces edge information when the kernel size is large.

Could you suggest to me any kernel that can satisfy the three conditions above, and is more robust than Gaussian in noise reduction, while maintaining edge information? I found a modified kernel, but it was very difficult to implement. Thank you so much.

There are many "Edge Preserving" filters in the image processing world.
2 very popular would be:

The Bilateral Filter.
Anisotropic Diffusion Filter.

Google search on each of the terms will give you plenty of data and code samples of each.

For any specific question about them, let me know, I'd be happy to assist.


Symbol Timing recovery for modulation producing ISI I am interested in understanding why the common timing recovery algorithms function for modulation schemes which produce ISI.

For example, suppose you are receiving at the output of a matched filter a raised cosine pulse and you are interested in timing recovery. For simplicity lets assume BPSK(though I would be interested in QAM). Besides the ideal sampling times there will be ISI.

Most of the common symbol timing recovery algorithms seem to use symmetry in the received signal, either balancing the measured value at the left and right sides of the ideal sampling instant or use zero crossings.

In the case of a sequence of raised cosine pulses, these features no longer exist as the signals overlap.

How do timing recovery algorithms work in this case?

In general the ISI is never severe enough for the symbol timing recovery algorithm not to work, but it does degrade its achievable performance on the $P_{be}$ vs. $\dfrac{E_b}{N_0}$ curve.  In other words, for the same error rate, you need better SNR to make up for the ISI.

ISI will degrade the S-curve of the Timing Error detector, causing it to flatten a bit and reduce the TED gain.  Gaussian noise will do the same.


Image Edge Preserving Smoothing I'm searching for a method that can smooth a 3D volume whilst preserving the edges in my volume.
I researched anisotropic diffusion filtering and bilateral filtering, but I'm having trouble to evaluate if they can be efficiently used on +30M voxels volumes.
So I need:

Relatively fast algorithm even with 512512200 volumes
Good edge preserving properties
Any help to understand those methods.

Thanks

If you can use the Bilateral Filter then you can use the Guided Filter.
The nice property of the Guided Filter is its low complexity.
There is a simple and efficient implementation with with linear complexity of the number of pixels.
You may also have a look at:

Fast Anisotropic Curvature Preserving Smoothing.
Simple Image Edge Preserving Filter.


What is a good example of an ergodic process? I'm trying to find simple examples of an ergodic process. What process comes to your mind as a good illustration of its properties?

A quick research (Wikipedia, another answer) mainly gives examples of non-ergodic processes. Also, I'm wondering which real world phenomena lend itself to be modeled as an ergodic process?

Just suppose I give you a series of numbers, and I tell you they were picked randomly. And you know I am not trying to deceive you. Numbers are: $3$, $1$, $4$, $1$,  $5$, $3$, $2$, $3$, $4$, $3$.

I now propose you to predict the next one, or at least, to be as close as possible. Which number would you pick?

[Think]

[Compute]

I bet most of the readers are likely to choose a number between $0$ and $6$. Because of the limited span.
Perhaps an integer. Who is likely to propose $\pi$ (even thinking of the first digits)?
Possibly  $2$, $3$, or $4$. Maybe even $3$.

Basically, you are assuming that I provided numbers with some unknown rule. And perhaps, you may think (or make the hypothesis) that the series of given numbers, if long enough,  can provide you with a good understanding of the rules that I have in mind.  If you do so, you hypothesize that  my mental process is ergodic:

  a process in which every sequence or sizable sample is equally representative of the whole (as in regard to a statistical parameter) (Merriam-Webster)

Here, there is no way to be sure that my series follows an ergodic process. 3432 is my card PIN, 3 a mistake (I intended 6, but I am clumsy), 4, 3, 1 and 5 are first digits of $\pi$ that I use quite often. My next "number" would have been C (in hexadecimal). I do not believe this process is ergodic. Each number comes out from  different laws. But honestly, I do not know. Maybe I am subject to some higher order forces that drive me under ergodicity rules.

So, ergodicity is a hypothesis of a sort of "simplicity" in the rules of a process. Like stationarity or sparsity.
 Cast a regular die with $6$ faces. Toss a normal coin. If nothing outside tries to influence the result (an invisible being that catches the die and shows some face of its choice), you are likely to produce an ergodic process.

Instead of being able to toss an infinite number of coins, with your   infinite number of thumbs, precisely at the same second, you toss one coin every second, and believe the final result is about the same.

The Brownian motion possesses ergodic properties too.


Amplification of a speech signal in the recorded audio I have recorded an audio.
I dont know how it happened that only one sided speech is recorded and the other speech is recorded with a very low sound. Is there any solution to amplify the other side signal.  I have tried to amplify with AVS software but its also ampliying the noise in the audio.
any help would be much appreciated.

I wouldn't recommend amplifying the low signal because that will increase noise significantly. As was pointed out above, it's better to downmix to mono.
Open your recording in a DAW (e.g. Audacity), delete the quiet channel, then copy the loud signal into the deleted channel. Voila!


Why use cross-spectral density to calculate frequency response? I am taking a class about system identification and currently learning about cross-spectral density. My textbook says that the frequency response, G, of a system can be determined as
$$G=\frac{S_{uy}}{S_{uu}}$$
where $S_{uu}$ is the Fourier transform of the autocorrelation of the input, and $S_{uy}$ is the Fourier transform of the cross-correlation of the input and the output (the cross-spectral density).

Now, in basic control systems classes you lean that $G=\frac{Y}{U}$.

The former method, using the cross-spectral densities, seems to be related to this more basic calculation, but I'm unsure why we bother using the cross-spectral densities at all. Is it simply to help eliminate noise from the signals or is there some other reason for it?

Thanks!

in the audio biz, we call this the "Two-channel FFT".  the cool thing about it is that you can measure the magnitude response of a room or something using music (that is decently broadbanded) as the test signal (and divide the output spectrum by the input spectrum).  you don't need to pollute peoples' ears with an impulse train or a maximum-length-sequence or pink noise or with a swept-frequency chirp.


When Is a Kalman Filter Different from a Moving Average? this thread asks when a discrete time Kalman filter is better/different from a simple moving average of the observations:

Why use a Kalman filter instead of keeping a running average?

there's no definitive answer. can someone give a definitive example where the kalman filter, ideally in simple 1D case, does something different (and better) than keeping a moving average, and state the conditions when the kalman filter would reduce to a simple moving average?

one thought is that the kalman filter would not weigh all the data points equally because its variance is initially smaller and gets better with time. but it sounds like that would only matter close to initial observations and that once the variance converged, the kalman filter would weigh each observation equally just like a moving average, so don't see when the two are different and why/when the filter would do better.

as the first answer (with the most votes) says, the kalman filter is better in Any case when signal is changing.
Notice the problem statement These use the algorithm to estimate some constant voltage. How could using a Kalman filter for this be better than just keeping a running average? Are these examples just oversimplified use cases of the filter? using a kalman filter to estimate a constant voltage is definitely, overkill. In that particular problem it is better to use the running average, which we know is the best estimator for Gaussian distributions. In this example the measured voltage is the actual voltage V but with some noise typically modeled as 0 mean Gaussian (white noise). so that our measurements are Gaussian with mean=V, and sigma=sigma noise.
The kalman filter is better suited for estimating things that change over time. The most tangible example is tracking moving objects. Let's imagine throwing a ball, we know it will make a parabolic arc, but what will our estimators show?
A Kalman filter will be very close to the actual trajectory because it says the most recent measurement is more important than the older ones (when the covariance is low that is). The running average takes all measurements equally

blue-ball trajectory, red-running average
(sorry no kalman if i have time i'll throw it in there if i have time, but it would me much closer to the blue line assuming you modeled the system well)

equations courtesy of wikipedia
The kalman filter on the other hand says, if our convariance and residual were small (meaning we had a good estimate), then we are going to stick with out previous estimate and tweak it a tiny bit based on the residual(or our estimation error). Now since our xhat k|k is very close to the actual state, when we make out next update, we will use a system state that closely matches the actual state.
At x=30, the running average says, the initial condition y(0) is just as important as y(29), thats that, and you get a huge error. The kalman filter accounted for this. It said since our error last time was huge, lets make a drastic change to our estimate (our xhat) so when we use it for the next update, it will be closer to what is actually happening
I hope that makes some sense
-------------------------EDIT--------------------------------
I just noticed your question asks about a moving average vs kalman. I answered running avg vs kalman (that is the topic of the link you provided)

Just to add a little more info specifically to the moving(windowed) average. The moving average is a better estimator of changing values. Since it only takes into account more recent samples. Unfortunately, it has a lag associated with it, especially around changing derivatives (Just look near t=30, where the derivative is going from positive to negative). This is because the average is slow to see fluctuation. Which is typically why we use it, to remove fluctuation (noise). The window size also plays a role. A smaller window is usually closer to the measured values, which makes sense and sounds good, right? The downside of this is if you have noisy measurements, a small window means more noise shows up more in the output. Let's look at the other question again

measurements with mean = .5, sigma = .1;
z = [0.3708435, 0.4985331, 0.4652121,...

the average of the first 3 samples is 0.4448629 not exactly close to the .5 expected value. This again shows, that with the smaller window, noise has a more profound effect on the output.
So then logically our next step is to take larger windows, to improve our noise immunity. Well, turns out larger windows are even slower to reflect actual changes (again look at t=30 in my graph) and the most extreme case of windowing is basically the running average (which we already know is bad for changing data)
Now back to the magical kalman filter. If you think about it it similar to a 2 sample windowed average (similar not the same). Look at X k|k in the update step, it takes the previous value, and adds to it a weighted version of the current sample. You might think, well what about noise? Why isn't it susceptible to the same problem as windowed average with a small sampling size? Because the kalman filter takes into account the uncertainty of each measurement. The weighting value K (kalman gain) can be though of as a ratio between the covariance(uncertainty) of your estimate and the covariance (uncertainty) of the current estimate (actually its the residual, but it's easier to think of it this way). So if the latest measurement has lots of uncertainty K decreases, and thus the most recent sample plays a smaller roll. If the latest measurement has less uncertainty than the prediction, k increases, and now the new information plays a bigger roll in the next estimate. So even with a small sample size, the kalman filter is still blocking out a lot of the noise.
anyway, I hope that answers the "windowed avg vs kalman" question now


ISI Nyquist criterion and spectral efficiency While reading Bernard Sklar's Digital Communications book I've got confused with ISI Nyquist criterion (not Nyquist rate criterion). There is stated (Sklar, 3.3 section) for symbol rate $R_s$ minimum bandwidth $W = R_s/2$ and for practical cases we have equation:

$W = 1/2 \cdot (R_s + r \cdot R_s)$,

where $r = 0...1$. So theoretically we have 2 symbols per Hz. For BPSK it is 2 bits/s per Hz. For QAM-64 it is 12 bits/s per Hz. So this statement I can't understand, because through my practice I used to think if we use BPSK and we want to achive e.g. 6 Mbps we need to allocate about 6 MHz of bandwidth (actually more because of rounding and synchronization/code redundancy). Where is my mistake in understanding of Nyquist theory?

Upd. Here is the picture from the book, where Nyquist channel is pictured

It seems the bandwidth of such a channel is calculated from $0$ to $1/2T$ and negative frequencies are neglected. Why in the baseband case we need to determine bandwidth in such a way? This can be an answer for the question above.

I don't have Sklar's text but I guess he is referring to the baseband case. For bandpass signals the required bandwidth is doubled. So in the bandpass case you get a maximum spectral efficiency of 1 bit/s/Hz for BPSK, and $\log M$ bits/s/Hz for $M$-ary QAM.

Just as a side note, spectral efficiency is measured in bits/s/Hz (and not in bits/Hz as you seem to suggest).


Phase Information at Higher Frequencies in Continuous Wavelet Transform I'm using the code I found here to compute the wavelet transform of a sine wave with a constant frequency.

#!/usr/bin/python2

from pylab import *
import matplotlib.pyplot as plt
import numpy as np
import scipy

x = np.linspace(0, 10, 65536)
y = np.sin(2 * pi * 60 * x)
N = len(y)
Y = np.fft.fft(y)

J = 128
scales = np.asarray([2 ** (i * 0.1) for i in range(J)])

morletft = np.zeros((J, N))
for i in range(J):
    morletft[i][:N/2] = sqrt(2 * pi * scales[i]) * exp(-(scales[i] * 2 * pi * scipy.array(range(N/2)) / N - 2) ** 2 / 2.0)

U = empty((J, N), dtype=complex128)
for i in range(J):
    U[i] = np.fft.ifft(Y * morletft[i])

plt.imshow(abs(U[:,scipy.arange(0,N,1)]), interpolation='none', aspect='auto')
plt.title("Sine Wave")
plt.xlabel("Translation")
plt.ylabel("Scale")
plt.show()

The result looks alright.

What I'm really interested in looking at is the phase. I'm extracting phase from the above code using:

imshow(np.unwrap(np.angle(U)), aspect='auto')

and it looks like this:

Why is the phase information present at frequencies higher than (or scales lower than, NOTE: inverted y-axis) that of the signal?

As you can see there is energy also in higher frequencies than your original.
This is probably since every level of the wavelet has a band, so the wavelet can show energy at frequencies higher.

The phase exists only as long as there is energy... so it directly connected to it. You can see if the data on the spectrum exists only where there is non zero energy in the spectrum.


What is Label Refinement in the Context of Image Segmentation I have been reading the paper entitled Measures of Similarity by Unnikrishnan and Herbet, in which they discuss some metrics to evaluate image segmentation. They use the word refinement of labels/segmentations/images in different places. However, it is unclear to me what this means, and an online search for "image refinement" does not bring much enlightenment.
I hope someone can elaborate a bit on the meaning of refinement in the context of images, or more specifically if possible, in the context of image segmentation.
Thank you in advance!

The correct context of the refinement key word is segmentation.
Label Refinement in the context of image segmentation is a step to increase the resolution and understanding of the segmentation.
It can be done by exterior knowledge (Like labels on features) or other optimization steps to have a better results of the segmentation (Which basically labeling of the image pixels).
For instance, if the segmentation is based on Super Pixels a Refinement Step could be merging of some Super Pixels.
You may have a look at Label Refinement Network for Coarse to Fine Semantic Segmentation. In this paper they use Multi Scale approach for label segmentation. By using coarse level to refine finder level.


Determine when a constant signal changes to another constant value As a novice to digital signal processing I have a scenario that probably seems trivial but I dont know the answer to. I want to detect when a signal changes but ignore noise.

  255 255 255 255 255 210 180 170 140 130 130 130 130 130 130

In the above signal there is a constant 255 value then it changes to a constant 130 value (with some noise in between). Is there a signal processing technique that can detect this, ie, it ignores noise and just detects the constant signal? Note the noise can persist for a while but we will always eventually get back to some constant signal value.

Context: I am learning computer vision as a hobby. My example signal represents grayscale intensity values (0-255). I want to detect the mode change, ie, detect when the most common value being seen changes whilst ignoring noise. I am aware of image gradients and gradient magnitudes identify areas of rapid signal change. However my signals wont always have a rapid change, they will always however have constant signal values between the change/noise. But this noise could be short or very long. Also the difference between the constant values could be very small or very larges. So it could go from constant value 255->240 or 255->0.

As @AlexR mentioned, the digital 1st order derivative is the most straight forward and probably the most trivial solution to the problem, depending on what your signal looks like and the rate of changes involved between samples. The derivative can then be processed in a sliding window to find where it transitions and where it settles (derivative of constant will be close to zero). With these techniques you can reliably detect constants ONLY if the derivative of the signal is greater than the noise amplitude before the signal settles to its constant value. If the rate of change of the signal is slow (i.e. the sample to sample transitions are on the order of magnitude of your noise) then it will be difficult to detect the constant signal from this method. I scale the signals below to be on the same order of magnitude and also take the absolute value of the derivative.

close all
clear all
%signal
x = [ones(1,5)*255 253 251 ones(1,6)*240 240 230 220 200 193 187 180 170 164 140 130 130 120 110 80 ones(1,10)*70 80 95 112 125 137 140 170 180];

%repeat signal
y = repmat(x, [1 5]);

%max value 2^8
MAX = 255;

%noise amplitude
Namp = 5;

%add noise to signal
y = y+randi(Namp,1,length(y));

%1st order derivative
c = diff(y,1);

%plots
figure
stem(abs((c)/MAX),'color','b')
hold on
plot(y/max(y))
line([1,length(y)],[0,0],'color','r','linewidth',2)
line([1,length(y)],[Namp/MAX,Namp/MAX],'color','r','linewidth',2)
ww = 4
reg = zeros(1,ww)
for i = 1:length(c)-1
  reg = [reg(2:end) abs(c(i)/MAX)];
  if(reg(1) >Namp/256 & sum(reg<Namp/256)>= 3 )
      line([i-ww+2, i-ww+2],[0 1],'color','g','linewidth',2)
  end
end
title({'Red: Signal', 'Blue: Derivative', 'Green: start of constant signal'})
hold off

The signals generated from the derivative algorithm can be used to make a clean piecewise continuous function as seen below. The delay of which, will be a function of your window size-1 sample.

Updated code:

close all
clear all
%signal
x = [ones(1,20)*255 253 251 ones(1,4)*240 255 255 255 255 ones(1,5)*240 240 230 220 200 193 200 220 ones(1,10)*230 210 190 187 180 170 164 140 130 130 120 110 80 ones(1,10)*70 80 95 112 125 137 140 170 180];
%x = [ones(1,5)*130 ones(1,10)*255 ones(1,10)*20];
%repeat signal
y = repmat(x, [1 5]);

%max value 2^8
MAX = 256;

%noise amplitude
Namp =5;

%add noise to signal
y = y+randi(Namp,1,length(y));

%1st order derivative
c = diff(y,1);
c2 = diff(y,4);

hold off
%plots
figure
stem(abs((c)/MAX),'color','b')
plot(abs(c2)/MAX,'color','r')
hold on
plot(y/max(y))
line([1,length(y)],[0,0],'color','r','linewidth',2)
line([1,length(y)],[Namp/MAX,Namp/MAX],'color','r','linewidth',2)

ww = 4;
reg = zeros(1,ww);
filtered = []
piecewise =y(1);
for i = 1:length(c)-1
  reg = [reg(2:end) abs(c(i)/MAX)];

  %if we detect a transition to 0 derivative
  if(reg(1) >Namp/256 & sum(reg<Namp/256)>= 3 )

      line([i-ww+2, i-ww+2],[0 1],'color','g','linewidth',2)
      piecewise = y(i);
  end
  %if we detect a transition out of 0 derivative
  if( sum(reg<Namp/256)>= 3 & reg(end) >Namp/256 )
      line([i, i],[0 1],'color','y','linewidth',2)
      %piecewise = y(i)
  end

  filtered = [filtered piecewise];

end
title({'Red: Signal', 'Stem Plot: Derivative', 'Green: start of constant signal', 'Yellow: end of constant signal'})
hold off

figure
plot(filtered(3:end))
hold on
plot(y)
title({'Red: Signal', 'Blue: Cleaned Signal'})

If we take $x[n]$ to be the input and L to be the window length. First we take the 1st order digital derivative of $x[n]$ to get $\hat{x}[n]$

$$\hat{x}[n] = x[n+1]-x[n]$$

and then we examine $\hat{x}[n]$ in an L length window, noting that when the oldest sample in the window is large (signal is transitioning) and the rest of the samples in the window have settled close to zero (signal is constant) we are at an area of interest (we have found a constant signal after a transition). The longer the window L, the more samples you are considering after the transition before you trigger your algorithm. That is, if you know the signal will settle to a constant for a long period of samples, then choose a big L. If you want the output to transition quickly though, choose a smaller L. We know the smallest transition is zero but will be on the order of $N_{amp}$, the amplitude of the noise. If we scale everthing to the maximum amplitude the signal can be, $A_{max}$, then the newest L-1 samples are on the order of $\frac{N_{amp}}{A_{max}}$ and the oldest sample $x[n-L]$ is some value above that by $\epsilon$

We need to capture this mathematically to say: when the first $L-1$ samples are less than or equal to the noise amplitude and the sample L before is larger than the noise amplitude we have gone from a transition to a constant and the output should be the sample at this instance. Intuitively it seems that:

$$\frac{1}{(L-1)} \sum_{i=0}^{L-1}\frac{|\hat{x}[n]|}{A_{max}} \leq \frac{N_{amp}}{A_{max}}$$

should capture some essence of part of this, but doesn't quite capture that ALL samples in this $L-1$ window should be less than or equal to the noise amplitude.

it's tricky to capture this if-else statement mathematically. One way to get close to a single equation is to say the areas of interest (green lines in my plot) occur when

$$r[n] = \frac{\hat{x}[n-L]}{\frac{1}{(L-1)} \sum_{i=0}^{L-1}|\hat{x}[n]|} > 1 + \frac{\epsilon}{N_{amp}}$$

since the division above expresses:

$$\frac{\frac{N_{amp}+\epsilon}{A_{max}}}{\frac{N_{amp}}{A_{max}}}$$

This works ok, for example if you say you are only interested in $\epsilon$ greater than 5x your noise amplitude, then you get close, but it misses some transitions. It also seems to behave oddly when the window increases.

That is to say

$$y[n] = \left\{
\begin{array}{ll}
  x[n], \quad & r[n] > 1 + \frac{\epsilon}{N_{amp}} \\
  y[n-1], \quad & otherwise \\
\end{array}
\right. $$

code to reflect this math:

    close all
clear all
%signal
xx = [ones(1,50)*255 253 251 ones(1,50)*240 ones(1,4)*255 ones(1,10)*240 240 230 220 200 193 200 220 ones(1,20)*230 210 190 187 180 170 164 140 130 120 110 80 ones(1,25)*70 ones(1,2)*256 ones(1,25)*70 80 95 112 125 137 140 170 180];
%x = [ones(1,5)*130 ones(1,10)*255 ones(1,10)*20];
%repeat signal
x = repmat(xx, [1 5]);

%max value 2^8
MAX = 256;

%noise amplitude
Namp =5;

%add noise to signal
x = x+randi(Namp,1,length(x));

%1st order derivative
xhat = diff(x,1);

L = 5;
reg = zeros(1,L);
filtered = [];
y = [];
y(1:L) =x(1:L)

for n = L+1:length(xhat)-1
  LminusOneSum = 0;
  for i= 1:L-1
    LminusOneSum = LminusOneSum + abs(xhat(n-i))/MAX;
  end
  r(n) =  abs((xhat(n-L))/MAX)/(LminusOneSum*(1/(L)));

  if(r(n) > (1+6))
      line([n-L+2, n-L+2],[0 1],'color','g','linewidth',2)
      y(n) = x(n);
  else
      y(n) = y(n-1);
  end

end
title({'Red: Signal', 'Stem Plot: Derivative', 'Green: start of constant signal', 'Yellow: end of constant signal'})
hold off

figure
plot(y(L:end))
hold on
plot(x)
%plot(r)
title({'Red: Signal', 'Blue: Cleaned Signal'})


can measurements go out of the uncertainty bounds? In the below picture, the measurements are inside the $\pm 3 \sigma$ bounds. In my experiment, the measurements sometimes go out of the uncertainty bounds.

This is a snapshot of my plot where the red is $\pm 3 \sigma$ bounds (i.e. $true \pm 3\sigma$ ), the black is the measurements, the blue is the true value and the green is the EKF. Is there something wrong with my results?

The meaning of $±3σ$ is that for a given distribution of the measurements, most of the measurements (which have values of up to three time the standard deviation around the average) will be in that boundary. By definition there will be values out side the boundary.

You need to make some specific computations to find out if in your implementation there is a deviation from your expected performance.

I wonder what is the distribution of the values you assumed? How you determine for each measurement the average and the boundaries - you have a kind of dynamic behavior which places a challenge on defining the various values.


Clarifying matrix notation from an ICA-CMN paper In the paper I am referring (and here from citeseer), complex vectors $\mathbf{z}$ and matrix $\mathbf{M}$ were defined as follows

\begin{align}
{{\bf z}} &= \left[z_{1},z_{2},\ldots,z_{N}\right]^{T} \in {\mathbb{C}}^{N} \\
{\tilde{\bf z}} &= \left[z_{1},z^{*}_{1},\ldots, z_{N},z^{*}_{N}\right]^{T} \in {\mathbb{C}}^{2N} \\
{\hat{\bf z}} &= \left[z^{R}_{1},z^{I}_{1},\ldots, z^{R}_{N},z^{I}_{N}\right]^{T} \in {\mathbb{R}}^{2N}\\
{\bf M} &\in {\mathbb{C}}^{N \times N}\cr \tilde{{\bf M}} & \in {\mathbb{C}}^{2N \times 2N} \cr {\hat{\bf M}} & \in {\mathbb{R}}^{2N \times 2N}
\end{align}

However, I can't seem to figure out how $\tilde{\mathbf{M}}$ would look like.
I have\begin{equation}
\mathbf{M} =
    \begin{bmatrix}
    z_{11} & z_{12}   & \dots     & z_{1N}\\
    z_{21} & z_{22}   & \dots     & z_{2N}\\
    \vdots & \vdots     & \ddots    & \vdots\\
    z_{N1} & z_{N2} & \dots         & z_{NN}
    \end{bmatrix}
\end{equation}
and the paper later clarifies
\begin{equation}
\hat{\mathbf{M}} =
    \begin{bmatrix}
    z_{11}^R & -z_{11}^I & \dots    & z_{1N}^R & - z_{1N}^I\\
    z_{11}^I & z_{11}^R & \dots    & z_{1N}^I &  z_{1N}^R\\
    \vdots   & \vdots    & \ddots   & \vdots & \vdots\\
    z_{N1}^R & -z_{N1}^I & \dots    & z_{NN}^R & - z_{NN}^I\\
    z_{N1}^I & z_{N1}^R & \dots    & z_{NN}^I &  z_{NN}^R\\
    \end{bmatrix}
\end{equation}
but not for  $\tilde{\mathbf{M}}$.
Thank you very much!

The matrices $\hat{\mathbf{M}}$ and $\tilde{\mathbf{M}}$ are constructed in such a way that the relation $\mathbf{M}\mathbf{x}=\mathbf{y}$ implies $\hat{\mathbf{M}}\hat{\mathbf{x}}=\hat{\mathbf{y}}$ and $\tilde{\mathbf{M}}\tilde{\mathbf{x}}=\tilde{\mathbf{y}}$.

Consequently, for constructing the matrix $\tilde{\mathbf{M}}$,  each element $m_{kl}$ of $\mathbf{M}$ must be replaced by a $2\times 2$ sub-matrix

$$\begin{bmatrix}m_{kl}&0\\0&m_{kl}^*\end{bmatrix}$$

In a similar way, the matrix $\hat{\mathbf{M}}$ is formed by replacing each element $m_{kl}$ of $\mathbf{M}$ by the $2\times 2$ sub-matrix

$$\begin{bmatrix}m_{kl}^R&-m_{kl}^I\\m_{kl}^I&m_{kl}^R\end{bmatrix}$$

where the superscripts $^R$ and $^I$ denote real and imaginary parts, respectively.


Conversion from counts to dB I have a sound signal, recorded by a hydrophone, which units are counts. I would like to convert the counts into dB re uPa (Decibels relative to microPascals), is that possible?

As a lot of people said in the comments if you have the sensitivity of your hydrophone (which you state that you do) and consider you are using your equipment in its linear regime, then you should be able to find the value in $\mu Pa$. From that, you could easily convert to $dB$.

The process is quite simple. I will illustrate an example here, taken from a real-life sound level survey equipment. It is quite normal for such microphones to exhibit a sensitivity of the order of $30 mV/Pa$. This means that when you measure the voltage output you can convert directly from voltage to pressure. As an example, let's assume that we measure at the output of the microphone an RMS voltage value of $108 mV$. From the sensitivity of the microphone, we know that for a pressure of $1 Pa$, the output should be $30 mV$, so now we can divide our output by the microphone's sensitivity (if you are not sure how we end up with this result you can look for the rule of three, which is exactly how you reach that) to get the pressure in $Pa$ that corresponds to the measured voltage. That is

$$Pressure = \frac{V_{out}}{sensitivity} \implies
Pressure = \frac{108 mV}{30 \frac{mV}{Pa}} \implies
Pressure = 3.6 Pa$$

Now, keep in mind that here, you have to take into account the possible gain of a preamplifier. That is if there is equipment with gain other than unity after the microphone you have to compensate for that. For example, in our case, if there was a pre-amplification stage after the microphone with a gain of 2, then we should have used a value of $\frac{108 mV}{2} = 59 mV$, which corresponds to the output of the microphone. This could also be done by subtracting the gain in $dB$ from the result (see next step).

Next step from that is to use the appropriate reference to convert to an absolute decibel value. In the air, the most widely used reference value is $20 \mu Pa$ (which is considered to be the hearing threshold for a healthy average person for a steady tone with a frequency of $1 KHz$). So, using that and taking into account that pressure is a derived, field quantity we get (SPL here stands for Sound Pressure Level)

$$SPL = 20 \log_{10} \left( \frac{p}{p_{ref}} \right) \implies
SPL = 20 \log_{10} \left( \frac{3.6 Pa}{20 \cdot 10^{-6} Pa} \right) \implies \\
SPL = 20 \log_{10} \left( 0.18 \cdot 10^{6} \right) \implies
SPL \approx 20 \cdot 5.255 \implies
SPL \approx 105.1 dB ~ re ~ 20 \mu Pa$$

Again, if we had an element with gain other than unity, this could be subtracted from this result. For example, a gain of 2 like above corresponds to a gain of $6 dB$ which we should subtract from this result (if not done in the previous step). This would give a result of about $105.1 - 6 = 99.1 dB_{SPL}$.

In the exact same way, you can calculate your quantities but instead of $mV$, you will use counts (since this is the quantity you have and the sensitivity is given in this quantity too).

Now, if you have the sensitivity in $dB$ as you state, all you have to do is to first convert that to $counts/\mu Pa$ and use that in the same way we used $mV/Pa$ sensitivity. In order to convert from $dB$ to linear scale, you use the properties of the logarithm and solve like

$$dB = 10 \log_{10} \left( x \right) \implies
\frac{dB}{10} = \log_{10} \left(x \right) \implies
x = 10^{\frac{dB}{10}}$$

where x is the value you are seeking (the linear value corresponding to the $dB$ value shown in your sensitivity rating), $dB$ is the sensitivity value you have. Now, you have to keep in mind that instead of the multiplier of the logarithm being $10$, it might be 20 (I haven't encountered any other value so far, but I don't know the counts as a unit either, so I can't say for sure) depending on the way it is derived.

Hope this helped somehow.


Band Limited Table Look Up Oscillator C++ I would like to create a Band Limited Table Lookup Oscillator that functions much like the sin() function in C++ works where i would inout a phase value and receive a sample value. for example Sample s = saw(x); or Sample s = square(x);.

To fill the table I would first generate an oversampled table by a factor of "X" that contains the wave in a non band limited state. Then I would run the table through a filter of some sort FIR or IIR most likely. Then i would decimate the table by a factor of "X" and hopefully then have a table that i could read from that would produce a band limited signal. I would also have to interpolate across the table to ensure that all phase angles will be accommodated for.

I would like to know if this theory would work? Meaning would this produce a signal that is Band limited.

I will mock up the basic code below.

double saw(double phase){
    return interpolate(phase);//i know this is incomplete
}

double fillTable(){
    double table[size*oversampling];

    //fill table here. skipping for length of example

    for(int i =0;i<tableSize;++i){
        table[i]=filter(table[i]);
    }
    double table2[size];
    int x=0;
    for(int i =0;i<tableSize;++i){
        table2[i]=table[x];
        x+=oversampling;
    }
}

//usage
sample x = saw(phase);

I know the code is incomplete and syntactically incorrect but i just wanted to give an example of what i would like to be able to do.

So again. Will this be Band limited?

Please note that you cannot implement a band-limited oscillator with a pure (stateless) function like:

Sample s = saw(x);

The reason is that the function must be aware of the rate at which x change to generate a signal with the correct bandwidth

For example, let us assume your sample rate is 48kHz. If you call your function in this order:

saw(0.00000000);
saw(0.00000001);
saw(0.00000002);

You can see that the argument is slowly changing to the point that aliasing is not an issue, so your function could return something like x / M_PI - 1 (assuming input is in [0, 2 pi[ ; and output in [-1, 1]) without any lookup table. But now, if you call it this way:

saw(0.0);
saw(3.0);
saw(6.0);

We are getting close to Nyquist frequency and the output must then be heavily band-limited - actually closer to sin.

So a band-limited sawtooth or square cannot be generated by mere waveshaping of the naive saw. You need your function to be aware of the target frequency so that it can make the right decision.

Another related misconception that can be found in your code is that you assume that you would need only one table. There is just not one single "band-limited sawtooth" waveform that would work for all frequencies. For low frequencies, there is very little different between the band-limited saw and the naive saw. For high frequencies, the band-limited saw is effectively a sine. So you have to define how much memory you want to dedicate to your task, and generate many band-limited tables corresponding to various frequencies. At run-time, you decide which one to use depending on the frequency of the signal. If this is for a musical application in which you don't want to hear the transition from one table to the next as notes rise in pitch, you'll need to add one further layer of interpolation on top of that to crossfade from one table to the next. A typical implementation would use one table by octave or per 8 semitones.


What is phase shift? I'm reading about when a $H(z)$ system is min phase or max phase, and I get when a system is a min phase or max phase depending on the location of the zeroes. But I don't get at all the phase shift part. on one point wikipedia says that these two transfer functions:
$\frac{s+10}{s+5}  \; and \;\; \frac{s-10}{s+5} $
$\\\\$
"have equivalent magnitude responses but however the second system has a much larger contribution to the phase shift." So what is phase shift?

Phase shift is often confused with time delay, although related these are two different quantities. Phase shift specifically refers to the rotation of a complex number. For example given the complex number $x = Ae^{j\theta}$ which through Euler's formula has real and imaginary components given as $A\cos(\theta)+jA\sin(\theta) = I + jQ$, a phase shift of $\phi$ would rotate $x$ according to:
$$ x e^{j\phi}= Ae^{j(\theta+\phi)}$$
Resulting in real and imaginary components as:
$$xe^{j\phi}= A\cos(\theta+\phi)+jA\sin(\theta+\phi)$$

Time domain and frequency domain waveforms can be represented as complex numbers where each sample has a magnitude and a phase. A phase shift would be a rotation of those samples, in either domain (so we can have a phase shift in time, or a phase shift in frequency). A constant delay in time has the property of a linear phase in frequency: the phase shift versus frequency will be linearly increasing in the negative direction for a fixed delay in time. For instance, given the observation of a sinusoidal tone at the input and output of a cable (fixed delay), for an extremely low frequency we would see the result of a very small phase difference between the sinusoid measured at the output compared to the input, but as the frequency is increased, this measured phase would increasingly larger in direct proportion to the increase in frequency.
For further related details, see:
What is an Intuitive Explanation of the Phase of a Signal
How does a shift in time domain result in phase shift in frequency spectrum?
Meaning of Phase response of a filter? In simple words?


What is the advantage of using Bessel function as the basis in a Hankel Transform? I came across literature which says that the Fourier - Bessel (Hankel) transformation is better than a STFT or a Wigner - Ville distribution as the Bessel function forms an orthonormal basis and they occupy finite bandwidth. This is supposed to translate into better time and frequency resolutions. Isn't it the same with the regular fourier transform ? so why not use an STFT itself to transform the signal into time frequency domain ?

The Hankel Transform includes the Bessel function as part of its definition. I wouldn’t say the Bessel has an advantage. I would say that the Hankel Transform has an advantage where a two dimensional signal has radial symmetry, or saying it another way, the signal in $x$ and $y$ can be expressed as a product of two functions in $r$ and  $\theta$ through a change of variables.
$x=r\cos \theta$ $y=r\sin \theta$

$$f(x,y)dx dy \Rightarrow g(r) h( \theta) r dr d \theta $$

The Hankel Transform is equivalent to a Fourier Transform on the radius $r$ variable.

Any advantage over STFT or Wigner would need to be explained in the literature you read.  Im not aware of any general use of Hankel grams, or why the Hankel would make any obvious relevance outside of two dimensional geometries with rotational symmetry.


Does DFT produces the same output as FFT? In my journey about learning what / why / how of DFT, I tried to implement a DFT on MATLAB and then I compared its output with fft output and then I noticed it was different, does DFT produces a different output from fft?

% it creates a simple gray image (4x4)
I = [255, 255, 30, 100
    255, 50, 90, 20
    70, 70, 20, 10
    100, 20, 10, 0];
% it converts it to grayscale
I = mat2gray(I);
% shows it
imshow(I);

IDFT = [];

% DFT
N1 = 4;
N2 = 4;
e = exp(1);

for w1=1:(N1)
  for w2=1:(N2)
      % X(w1, w2)
      DFT = 0;
        for n1=1:(N1)
            for n2=1:(N2)
                Intensity = I(n1,n2);
                FirstDimension = e ^ (-1i*( (2*pi/N1) * w1 * n1));
                SecondDimension = e ^ (-1i*( (2*pi/N2) * w2 * n2));
                DFT = DFT + (Intensity * FirstDimension * SecondDimension);
            end
        end
       IDFT(w1, w2) = DFT;
  end
end

IDFT
fft(I)


How can I identify Arabic letters in an image using Java? I want to identify Egyptian car number plates that look like this:

How can I do this in Java ?

I plan to do this in Android, as the user will catch the number plate with a camera and the results will be submitted to server.

You have to do following methods:

First step: Number plate extraction

You can use OpenCV for this purpose. It has Java api. Also it is ported to Android.

Second step : OCR

OCR is used to recognise the letters on notice board.

You can train your own OCR for arabic letters and digits using OpenCV.

Or you can use tesseract-OCR. I think it supports arabic also.

Following paper will be useful for you: Paper 1

Also, this site includes complete code for number plate extraction, but in C++. But it will be a good start.


Image Processing - find the centre of the estimated circle of eye I have been recently working on detecting the direction of eye pointed out and done some basic operations on the image of the eye. the below is the image of boundry of the eye.

Now the problem is how do i detect the centre of the rough circle i have obtained. please help me out by providing me a matlab code.

Unfortunately, I dont have access to a MATLAB image processing toolbox. I've coded a solution in python, showing the idea of how to solve it:

import skimage
import skimage.io
import scipy.ndimage as ndi
import skimage.morphology as m
import skimage.measure as measure

plt.figure(figsize=(8, 12))
I = skimage.io.imread("https://i.stack.imgur.com/eCVzl.jpg")
plt.imshow(I[:,:,2])
bw = I[50:250,100:300,1] < 100  # Import and binary image. Cut the region of interest (e.g. remove the white border)

plt.subplot(321)
plt.imshow(bw, interpolation='none')

plt.subplot(322)
plt.title('Distance transform of input image')
R, D = m.medial_axis(bw, return_distance=True)  # Calculate the distance transform (Matlab: bwdist)
plt.imshow(D, interpolation='none')

plt.subplot(323)
plt.title('Convex hull of input image')
C = m.convex_hull_image(bw == 0)  # Calculate the convex hull of the pixels (i.e. the are where the maxima can be)
plt.imshow(C)

masked = C * D  # zero-out all distances that are not within the convex hull
plt.subplot(324)
plt.imshow(masked)
plt.title('Masked distances according to convex hull')

M = measure.moments(C.astype(np.uint8), order=2);
cy, cx = M[0, 1] / M[0, 0], M[1, 0] / M[0, 0]   # Calculate the convex hull center of mass
x, y = np.meshgrid(np.arange(C.shape[0])-cx, np.arange(C.shape[1])-cy);
weights = np.exp(-0.003*(x*x+y*y))  # calculate a weighting mask

plt.subplot(325)
plt.imshow(weights)
plt.title('Distance weights')

masked_weighted = masked * weights

plt.subplot(326)
max_pos = np.unravel_index(masked_weighted.argmax(), masked.shape) # detect maximum pixel and its 2d-coordinates
plt.imshow(masked_weighted)
plt.plot(max_pos[1], max_pos[0], 'kx', mew=2)
plt.title('Weigthed distances with detected center')

plt.subplot(321)
plt.plot(max_pos[1], max_pos[0], 'kx', mew=2)
plt.title('input image with detected center')

plt.tight_layout()

The idea is based on an operation called "Distance Transform" and that the center of a circle is the point that has the largest distance to its boundary. The program assumes that the center of the circle is roughly in the middle of area spanned by the edges. Furthermore, it assumes that the maximum distance pixel corresponds to the circle center (i.e. other edges are closer together).  So, in rough steps:

calculate distance transform of Image
calculate convex hull of the edge pixels
caluclate center of the convex hull to detect the rough center
weight the obtained distances from 1) by the distance from the estimated center
use the strongest pixel as the center of the circle.

The code sohuld be easy to translate into Matlab, esp. as I have plotted all intermediate images.


Do $|s(t)|$ and $|S(f)|$ uniquely determine $s(t)$? Consider a signal $s(t)$. My question is if you know $|s(t)|$ and $|\mathcal{FT}[s(t)](f)| = |S(f)|$ or equivalently $|s(t)|^2$ and $|S(f)|^2$ is it possible to determine $s(t)$? That is, is $s(t)$ uniquely determined by $|s(t)|$ and $|S(f)|$?

My strong feeling is that the answer is no but I can't come up with two signals $s_1(t) \neq s_2(t)$ which have the same magnitude and whose Fourier transforms have the same magnitude. Actually I can come up with trivial examples but that doesn't fully satisfy me. Let me share by work so far.

First of all, if we consider

$s_2(t) = s_1(t) e^{i\phi}$

Where $\phi$ is any real constant number then we clearly see that

\begin{align}
S_2(f) = S_1(f) e^{i\phi}
\end{align}

and

\begin{align}
|s_2(t)| &= |s_1(t)|\\
|S_2(f)| &= |S_1(f)|
\end{align}

So this example is a clear counterexample. However, this is just a trivial global phase that cannot be determined. That is, it is clear now that there is global phase multiplier that cannot be determined if only $|s(t)|$ and $|S(f)|$ are known. But, lets say that we either don't care about an overall phase factor or we happen to know that $s(0) = |s(0)|$ so that this global phase can be fixed.

Is it possible to come up with a non-trivial example of two functions which are a counterexample to what I am asking? I can phrase it this way.

If $|s_2(t)|= |s_1(t)|$ then

\begin{align}
s_2(t) = \Phi(t)s_1(t)
\end{align}

Where $|\Phi(t)|=1$ for all $t$. We then have that (by the convolution theorem)

\begin{align}
S_2(f) = S_1(f) \ast \Phi(f)
\end{align}

The question is is it possible for

\begin{align}
|S_2(f)| = |S_1(f)\ast\Phi(f)| = |S_1(f)|
\end{align}

If a non-trivial $\Phi(t)$ could be demonstrated then it would be clear that $|s(t)|$ and $|S(f)|$ do not uniquely determine $s(t)$ in any sense but I can't seem to come up with anything.

If I am wrong and it turns out that $|s(t)|$ and $|S(f)|$ do in fact uniquely determine $s(t)$ then could anyone provide a method to determine $s(t)$ (up to a missing overall phase) from $|s(t)|$ and $|S(f)|$?

As pointed out in a comment by Jazzmaniac, the proof below is flawed. The mistake I made was the assumption leading to Eq. $(3)$. Of course, $|S_1(f)|=|S_2(f)|$ implies that the two spectra are related by an all-pass transformation, but the phase of that all-pass can be dependent on the input signal, and hence need not be LTI.

One example of two signals satisfying $|s_1(t)|=|s_2(t)|$ and $|S_1(f)|=|S_2(f)|$ without satisfying $(1)$ is $s_1(t)=e^{j\phi(t)}$ and $s_2(t)=e^{-j\phi(-t)}$. Clearly, $|s_1(t)|=|s_2(t)|=1$. Since $s_2(t)=s_1^*(-t)$, the corresponding spectra satisfy $S_2(f)=S_1^*(f)$, and, consequently, $|S_1(f)|=|S_2(f)|$.

The dual case would be $S_1(f)=e^{j\Phi(f)}$ and $S_2(f)=S_1^*(-f)=e^{-j\Phi(-f)}$. The signals $s_1(t)$ and $s_2(t)$ then satisfy $s_2(t)=s^*_1(t)$ and $|s_1(t)|=|s_2(t)|$.

The second example above requires the signals to be complex-valued (otherwise $s_1(t)=s_2(t)$ would hold). However, the first example also includes the real-valued case, i.e., a signal $s_1(t)$ that only attains two values: either $1$ or $-1$ (like a telegraph process). With $s_2(t)=s_1(-t)$ we guarantee $|s_1(t)|=|s_2(t)|$, and since $S_2(f)=S_1(-f)=S_1^*(f)$ we also have $|S_1(f)|=|S_2(f)|$.

OLD answer:

I think I can show that the trivial case

$$s_2(t)=e^{j\alpha}s_1(t),\qquad\alpha\in\mathbb{R}\tag{1}$$

is the only case for which $|s_1(t)|=|s_2(t)|$ and $|S_1(f)|=|S_2(f)|$ holds.

From the requirement $|s_1(t)|=|s_2(t)|$ it follows that

$$s_2(t)=e^{j\phi(t)}s_1(t)\tag{2}$$

with a real-valued function $\phi(t)$. From the requirement $|S_1(f)|=|S_2(f)|$ it follows that $s_2(t)$ can be obtained from $s_1(t)$ by all-pass filtering:

$$s_2(t)=\int_{-\infty}^{\infty}s_1(\tau)h_{AP}(t-\tau)d\tau\tag{3}$$

with $H_{AP}(f)=\mathcal{F}\{h_{AP}(t)\}$ satisfying $|H_{AP}(f)|=1$.

Since the multiplication in $(2)$ is memoryless, Equations $(2)$ and $(3)$ can only give the same result if $h_{AP}(t)=e^{j\alpha}\delta(t)$, i.e., for a trivial all-pass which just scales the input by a complex constant with magnitude equal to $1$. Consequently, the function $\phi(t)$ in $(2)$ must be constant (and equal to $\alpha$). This means that $(2)$ and $(3)$ are reduced to Eq. $(1)$.

Consequently, given $|s_1(t)|$ and $|S_1(f)|$, the signal $s_1(t)$ is determined up to a complex factor with magnitude $1$. I just do not know of any method to actually compute $s_1(t)$ from these two magnitude functions.


Shift estimation using fourier transform If i have two images of the same scene, reference image and a shifted image, the shift can be in x or y directions, i want to estimate the shift using the shift property of the fourier transform.

S1 = fftshift(fft2(im1)); % Fourier transform of the reference image
S2 = fftshift(fft2(im2)); % Fourier transform of the image to be registered
Q = S1./S2;
A = angle(Q);

How can i use the phase shift Ato estimate the shift in both directions?

Assume the image is 256x256, The shift in spatial domain corresponds to a linear phase in the frequency domain.

u=-128:127;
v=-128:127;
u=repmat(u,size(u,2),1);
v=repmat(v',1,size(v,2));

linear_phase= exp(-2*pi*1i.*(((v).*dy)+((u).*dx)));

dx shift in $x$ and dy shift in $y$, and A equals linear_phaseso i can get a lot of equations to solve for dx and dy

If dx or dy is $0.5/128$ the shift is one pixel in $x$ or $y$


Simplest way to generate AR(2) process on MATLAB As part of a project I need to use autocorrelation method of estimating model paramters of an autoregressive process on MATLAB.

Can anyone tell me the simplest way to generate an AR(2) process on MATLAB so that I can estimate its model parameters $\hat{a}_{p}(k)$ and $\hat{b}_{0}$?

The difference equation for the process is given below:

$$x(n) = -0.9x(n-1) + w(n)$$

The simplest way to approximate an AR-2 process in Matlab / Octave is the following:

N = 1024;                    % number of process samples.
a = [1, -0.9, 0.2];          % denominator coefficients, p = 2.
b = [1.0];                   % numerator coefficient.
x = filter(b,a, randn(1,N)); % generate N sample of AR-2 x[n].

Note: an AR process requires a true-white noise sequence $v[n]$ at the input of the filter but here we input a single instance of a crude approximation of it. Hence the process is not truly an AR-2 but an approximation...


why is the plot of the PSD of this autocorrelation function oscillating? I was doing an exercise trying to compute the PSD of the following discrete autocorrelation function:

\begin{equation}
r_x(m) =(1/2)^{|m|}
\end{equation}

In this simple  matlab code that I came up with to compute the PSD, I'm getting something that I don't quite understand, why is my frequency response oscillating?.If I plot the abs(Sx) i get a smooth plot but , the dft of a real and even sequence should be real and even as well in my understanding . Why I'm not getting a real frequency response and instead this oscillating response appears? Also, matlab is ignoring the imaginary part of the PSD computed as it says in the command window. What is going on  here?

N=100;
m=-N:1:N;
x=(0.5).^abs(m);
subplot(2,1,1)
stem(m,x)
Nfft=512;
Sx=fft(x,Nfft);
subplot(2,1,2)
plot((0:1:Nfft-1)/Nfft,Sx);

Here are  the figures:

You are right that the DFT of a real and even sequence is real and even. The problem is that your sequence is not even, at least not the way you programmed it. The fft command doesn't know how you defined the index m. By definition, the first data element in the vector x is associated with the index $n=0$:

$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N}$$

(where $N$ is the DFT length).

Your definition of the sequence x results in a complex spectrum, and the plot command plots the real part of it.

There are two solutions, one of which you already found: either you take the absolute value of the spectrum, or you define the data in such a way that they are indeed even. Due to the implicit periodicity of the DFT, this means that you append the left half of your data (for negative time indicies) at the end of your data vector. Note that the zero padding must now be done in the middle of your data vector:

x = [x(101:201),zeros(1,Nfft-2*N-1),x(1:100)];
Sx = fft(x,Nfft);
max(abs(imag(Sx))            % 2.4980e-16

The resulting spectrum has a small imaginary part due to numerical inaccuracies, so you always need to use real(fft(x,Nfft)), but it's always good to check if the imaginary part is really as small as it should be (around $10^{-16}$ for double precision floating point format, as used in Matlab).


Best way of segmenting veins in leaves? I have done a lot of research and found out methods like adaptive thresholding , watershed etc that can be used of detecting veins in leaves . However thresholding isn't good as it introduces a lot of noise

All my images are gray image please could anyone suggest what approaches to adopt while considering this problem in urgent need of help

EDIT:My original image

After thresholding

As suggested by the answer i have tried the following edge detection

Canny

Too much noise and unwanted disturbances

Sobel

Roberts

EDIT: Tried one more operation i get the following result its better than what i tried with canny and adaptive What do you feel?

You're not looking for edges (=borders between extended areas of high and low gray value), you're looking for ridges (thin lines darker or brighter than their neighborhood), so edge filters might not be ideal: An edge filter will give you two flanks (one on each side of the line) and a low response in the middle of the line:

ADD: If've been asked to explain the difference between an edge detector and a ridge detector more clearly. I apologize in advance if this answer is getting very long.

An edge detector is (usually) a first derivative operator: If you imagine the input image as a 3D landscape, an edge detector measures the steepness of the slope at each point of that landscape:

If you want to detect the border of an extended bright or dark region, this is just fine. But for the veins in the OP's image it will give you just the same: the outlines left and right of each vein:

That also explains the "double line pattern" in the Canny edge detector results:

So, how do you detect these thin lines (i.e. ridges), then? The idea is that the pixel values can be (locally) approximated by a 2nd order polynomial, i.e. if the image function is $g$, then for small values of $x$ and $y$:

$g(x,y)\approx \frac{1}{2} x^2 \frac{\partial ^2g}{\partial x^2}+x y \frac{\partial ^2g}{\partial x\, \partial y}+\frac{1}{2} y^2 \frac{\partial ^2g}{\partial y\, ^2}+x \frac{\partial g}{\partial x}+y \frac{\partial g}{\partial y}+g(0,0)$

or, in matrix form:

$g(x,y)\approx \frac{1}{2} \left(
\begin{array}{c}
 x & y
\end{array}
\right).\left(
\begin{array}{cc}
 \frac{\partial ^2g}{\partial x^2} & \frac{\partial ^2g}{\partial x\, \partial y} \\
 \frac{\partial ^2g}{\partial x\, \partial y} & \frac{\partial ^2g}{\partial y\, ^2}
\end{array}
\right).\left(
\begin{array}{cc}
 x \\
 y
\end{array}
\right)+\left(
\begin{array}{cc}
 \frac{\partial g}{\partial x} & \frac{\partial g}{\partial y}
\end{array}
\right).\left(
\begin{array}{c}
 x \\
 y
\end{array}
\right)+g(0,0)$

The second order derivative matrix $\left(
\begin{array}{cc}
 \frac{\partial ^2g}{\partial x^2} & \frac{\partial ^2g}{\partial x\, \partial y} \\
 \frac{\partial ^2g}{\partial x\, \partial y} & \frac{\partial ^2g}{\partial y\, ^2}
\end{array}
\right)$ is called the "Hessian matrix". It describes the 2nd order structure we're interested in.

The 2nd order part of this function can be transformed into the sum of two parabolas $\lambda _1 x^2 + \lambda _2 y^2$ rotated by some angle, by decomposing the Hessian matrix above to a rotation times a diagonal matrix of it's eigenvalues (Matrix decomposition). We don't care about the rotation (we want to detect ridges in any orientation), so we're only interested in $\lambda _1$ and $\lambda _2$

What kind of shapes can this function approximation have? Actually, not that many:

To detect ridges, we want to find areas in the image that look like the last of the plots above, so we're looking for areas where the major eigenvalue of the Hessian is large (compared to the minor eigenvalue). The simplest way to detect that is just to calculate the major eigenvalue at each pixel - and that's what the ridge filter below does.

A ridge filter will probably give better results. I've tried Mathematica's built in RidgeFilter (which calculates the major eigenvalue of the Hessian matrix at each pixel) on your image:

As you can see, there's only a single peak for every thin dark line. Binarizing and skeletonizing yields:

After pruning the skeleton and removing small components (noise) from the image, I get this final skeleton:

Full Mathematica code:

ridges = RidgeFilter[ColorNegate@src];
skeleton = SkeletonTransform[Binarize[ridges, 0.007]];
DeleteSmallComponents[Pruning[skeleton, 50], 50]

ADD:

I'm not a Matlab expert, I don't know if it has a built in ridge filter, but I can show you how to implement it "by hand" (again, using Matematica). As I said, the ridge filter is the major eigenvalue of the Hessian matrix. I can calculate that eigenvalue symbolically in Mathematica:

$\text{eigenvalue}=\text{Last}\left[\text{Eigenvalues}\left[\left(
\begin{array}{cc}
 H_{\text{xx}} & H_{\text{xy}} \\
 H_{\text{xy}} & H_{\text{yy}}
\end{array}
\right)\right]\right]$

=> $\frac{1}{2} \left(H_{\text{xx}}+H_{\text{yy}}+\sqrt{H_{\text{xx}}^2+4 H_{\text{xy}}^2-2 H_{\text{xx}} H_{\text{yy}}+H_{\text{yy}}^2}\right)$

So what you have to do is calculate the second derivatives $H_{\text{xx}}$, $H_{\text{xy}}$, $H_{\text{yy}}$ (using a sobel or derivative of gaussian filter) and insert them into the expression above, and you've got your ridge filter.


Express Circular Convolution in Terms of Linear Convolution g[n], h[n] : finite length sequences of length 7
yC[n]: circular convolution
yL[n]: linear convolution

Express yC[n] in terms of yL[n].

If same two sequences of equal length are convolved using both linear(by tabular method) and by circular(matrix method) the output is of differing length and values.

How can I possibly relate circular convolution to linear convolution?

$h[n]_N$ is the periodic extension of $h[n]$.

$$ h[n]_N = h[n + mN] \quad \text{for } 0 \le n+mN < N \text{ and } m \in \mathbb{Z}$$

so, given an $n$ and knowing $N$, you have to choose the correct integer $m$ so that

$$ 0 \le n+mN < N $$
or
$$ -n \le mN < N-n $$
or
$$ -\frac{n}{N} \le m < \frac{N-n}{N} = 1 - \frac{n}{N} $$.

Or, for the term $h[n-k]_N$,

$$ h[n-k]_N = h[n -k + mN] \quad \text{for } 0 \le n-k+mN < N \text{ and } m \in \mathbb{Z}$$

and

$$ -\frac{n-k}{N} \le m < 1 - \frac{n-k}{N} $$

so, you have to split your summation into two summations.


Phase correction algorithm: minimisation phase error Thank you for your participation in my discussion in advance.
I am working on the implementation of a phase correction.
Currently, I have finished the algorithm. I think it is a standard algorithm everyone of you uses.
Actually, it is the standard algorithm to correct a phase; I think.

Does anyone have an idea of how can I minimize the difference between angles of adjacent samples/symbols?

Had someone already an experience with such simulation and probably had still the simulation script and can share it with me?

I commend you for using an intuitive algorithm. However, there are already established algorithms with far better performance. Phase recovery algorithms work by filtering the error signal down to zero. An example is the Costas loop structure shown in the figure below for QPSK.

Let's start with phase error detectors (Highlighted in yellow). The arctan detector you used in your example works but it's not typically used in practical systems (The division of imaginary/real may result in overflow if the real quantity is very small). Most systems use the Costas and maximum-likelihood error detectors. Taking BPSK for example. The Costas phase error detector is simply ($real\times imag$). Error detectors are characterized by an S-curve obtained by finding the expectation of the error detector output (For BPSK should be the mean of ($real\times imag$)) for different phase error inputs (phase errors from $-\pi$ to $\pi$ in steps of $\frac{pi}{N}$) as shown below.

Filtering is done by something called the loop filter that can "drive" the error down to zero. The proportional-plus-integrator shown below is widely used at it can drive both static phase error and frequency offset to zero. $K_{n}$ parameters (which depend on the slope of the S-curve) dictate the behavior of the loop.

Regarding examples. There are plenty of examples online such as this one. You could find many more by searching for Costas loop examples online. If interested, you could have a look at GNU Radio tutorials. Plenty of examples too.
And finally, I could suggest some improvement in your simulation. You should consider adding a small frequency offset as well as phase noise to fully represent a real-world scenario.
Good luck.