Find similar patterns in a signal Black rectangles on these charts show same signal, which is actually color sequence from a image, taken from different positions.

It is easy to see that all images contain same pattern. Is here some algorithm which can detect these similarities, ideally implemented in Python ?

What about cut signals in multiple windows and try computing cross corelation on them ?

More detail like PAK-9 asks, and example images with the pattern indicated, would help get better answers.  If you update your question with details and data, you're also more likely to get example implementations.

That said, this looks to be a candidate for a wavelet transform (like with pywavelets), which would allow for the stretching and scaling that you need to capture the variations of the pattern. Pre-computed wavelets can be browsed here, or you can roll your own (like the fast 2D CDF 9/7 lifting transform).

You suggested cross-correlation, which might be good for a first pass (see scikit-image template matching for ideas).  You might modify this to loop through windows to account for the stretching and scaling, but you'll want to normalize them if you intend to subsequently use a classifier.  Also, the relative width and spacing of the peaks seem conserved, which would lend itself to windowed supervised classification by many means, like SVM or Random Forests (scikit-learn has plenty of classifiers and examples for you).


Why do we need to make simulation not only on baseband but also on passband? What is the purpose of doing passband simulation, what are the signal processing  imperfections we can reproduce in passband?
In other words what are we missing if xe simulate only in base band?

It is IQ imbalance, time/frequency estimate,  matching filtering ... ?

So, we'll have to look at two things: Passband Signals and Passband Systems.

Let's start with the general notion that the consideration in domain A is equivalent to  consideration in domain B if there's an invertible mapping between A and B that covers all possible inputs and outputs of the consideration.

Signals

That mapping for the complex baseband and real passband domain for signals is well known; mathematically, it's the Hilbert transform, followed by multiplication with a complex exponential of negative center frequency of the passband.

In the other direction, it's the multiplication with a complex exponentionential of a frequency high enough so the complete baseband is at positive frequencies, and then removal of the imaginary part.

That immediately posts three requirements:

The Hilbert transform must be well-defined
There must be a center frequency for the passband signal
There must be a frequency higher than half the bandwidth of the baseband signal.

2 and 3 are equivalent, as far as I can see – you can only use baseband if your signal has been a passband signal before¹.

About the 1. requirement:

Tough one, I've been trying to generally state a domain of definition for the Hilbert transform, but that's nothing for a lazy Sunday noon. Concretely, and function $f(t)$ that fulfills

$$\left\lvert\int\limits_{-\infty}^{\infty} \frac{f(t)}{t-\tau}\right\rvert < \infty$$

is Hilbert-transformable, and you need to account for the fact that you can use all your Lebesgue intergral magic here (you'll really have to, that thing usually has a singularity, and that is an example often used to motivate Cauchy's principal value theorem for electrical engineers).

Now, if $f\in L^p,\, p>1$, then that's a given; but periodic (and hence, not-energy-)signals do work too; tempered distributions do, and a few other things.

Considering the fact that Parseval's theorem already tells us that, inherently, a signal that's limited in energy in frequency domain (and every band-limited signal has to be) needs to be in $L^2$, we can infer that there is none band-limited passband signal that you cannot equivalently observe in baseband.

Systems

This is the tricky part. We're looking for passband systems that we can't map back and forth from passband to baseband.

¹ – that's true for more than signals with closed interval support; compactness suffices, as far as I can tell, but your mixing function then is no longer a single complex exponential. For example, you could have an infinite real tone comb with frequencies $N\pi,\,N\in\mathbb N$, and your mapping to baseband could be Hilbert followed by multiplication with a complex tone comb $\sum_{m\in\mathbb N} e^{jmt}$, and would still end up with a proper baseband, and a way of converting that back to "strangely defined passband". But that's more about generalizing the notion of pass- and baseband than it goes to answer the question...


Pulse Shaping in Digital communication I am new here. I hope I can get help.

I want to ask a question regarding digital communication.

Does anyone know what are pulse shaping used for? Are they used for ISI mitigation or for something else? Also what is the receiver side of this operation is it matched filter?

As you've told the main goal of pulse shaping filter is to restrict signal bandwidth in the transmitter while in receiver you should place the matched filter to maximize SNR. In the practical modems pulse shaping filter in the transmitter and matched filter in the receiver are usually the same in terms of frequency response. The task is to choose filters to minimize ISI (not all TX-RX filters pairs satisfy this condition while if even maximizing SNR value). From the theory is known raised cosine filter satisfy ISI condition. You can read about this filter here. Let's place such a filter in TX and RX side:

$H_{system} = H_{Tx} * H_{Rx}$

We can see the total response of the system is no longer raised cosine. So the root raised cosine filter was proposed. The main idea is

$H_{rrc} = \sqrt{H_{rc}}$, so for our system

$H_{system} = H_{Tx} * H_{Rx} = H_{rrc} * H_{rrc} = (\sqrt{H_{rc}})^2 = H_{rc}$

By doing this we can minimize ISI and maximize SNR if the filter bandwidth is closely matched to signal's one. So it is very practical and widely used in digital modems (especially single carriers (SC) modems).


Comb-like peaks in frequency spectrum I am currently working on a project which involves measuring the sound produced by the onset of a vapour bubble in boiling water. All the experiments went perfectly well and I'm now trying to obtain a frequency spectrum from the signal.

The sound-pulses produced by the bubbles only last about 40 ms max, simply too short to obtain a decent frequency resolution. That is why I'm measuring a numerous amount of almost identical pulses over a period of 5 seconds. The signal is being sampled at a rate of 100 kHz, which results in 500.000 samples. The interval between the sound-pulses is 95 ms ish as can be seen in the figure of the signal below:

EDIT: Added an additional figure of the signal:

When I create the frequency spectrum of the entire signal, I get this comb-like frequency spectrum with an interval of about 10 Hz between the peaks (which corresponds to the 95 ms time interval, but I would then rather expect a peak at 10 Hz in the spectrum):

I can't really get much information from this spectrum. I do know for instance that there should be a peak at around 167 Hz by looking at the time-domain (there's a small sine wave visible with about 6 ms peak to peak). The envelope of the spectrum does resemble the spectrum of a single sound pulse as shown below (but way less accurate):

The matlab code I used for the fourier transform was:

y_fft = abs(fft(y-mean(y)));
y_fft = y_fft(1:Nsamps/2);
f = Fs*(0:Nsamps/2-1)/Nsamps;    %Prepare freq data for plot

Does anyknow how to get a usefull spectrum in this situation or please tell me what I'm doing wrong?

Thanks in advance,

Marc

Your signals look quite similar to vowel signals - in which a periodic pulse train is filtered by a resonant filter.

I suggest you to look into AR modelling to decompose your signal into a residual, and a set of AR coefficients. The spectrum of the residual will give you the "pulse train" component of your signal (a sharp peak every time a bubble is produced) ; while the frequency response of the all-pole filter associated with the set of AR coefficients will contain information about the sound of each bubble. The 167 Hz peak you are looking will be easy to spot on the later.


Downsampling, shifting, high pass and low pass filter commutativity strong textI have been reading "The Stationary Wavelet Transform and some Statistical Applications" by Nason and Silverman, and there is a claim in the their paper of which I cannot convince myself.

Defininitions:

$\mathcal{H}$: A low-pass filter.

The action of the low pass filter on a doubly infinite sequence $\{..., x_{-1}, x_0, x_1, x2, ... \}$ is defined $(\mathcal{H}x)_k = \sum_n h_{n-k} x_n$
The filter satisfied internal orthogonality, i.e. $\sum_n h_{n} h_{n + 2j} = 0$

$\mathcal{G}$: A high-pass filter, defined by the filter $g_n = (-1)^n h_{1-n}$

Clearly $\mathcal{G}$ satisfied the same internal orthogonality, and is mutually orthogonal to $\mathcal{H}$, i.e. $\sum_n h_{n} g_{n + 2j} = 0$ $\forall j$

$\mathcal{S}$: A shift operator defined by $(\mathcal{S} x)_j = x_{j+1}$
$\mathcal{D}_0$:  A binary decimation operator that chooses every even member of a sequence.

$(\mathcal{D}_0 x)_j = x_{2j}$
$\mathcal{R}_0$ the inverse operation.

$\mathcal{D}_1$: A binary decimation operator that chooses every odd member of a sequence, similar to $\mathcal{D}_0$ above.
It is not even necessary for the same choice of "odd" or "even" to be used
throughout.

Suppose $\epsilon_{J-1}, \epsilon_{J-2}, ..., \epsilon_{0}$ are a sequence of 0's and 1's.
One can use the binary operator $\mathcal{D}_{\epsilon_j}$ at level $j$, and the original signal can be recovered by applying the corresponding sequence $\mathcal{R}_{\epsilon_j}$

On page 4, the paper claims,

  "If $x$ is a finite sequence, define the shift periodically at the boundary. It then immediate from the definitions that $\mathcal{D}_1$ = $\mathcal{D}_0 \mathcal{S}$ and hence that $\mathcal{R}_1$ = $\mathcal{S}^{-1} \mathcal{R}_0$.

NOW comes the portion I can't seem to grasp...

  It is also easy to see that $\mathcal{S} \mathcal{D}_0 = \mathcal{D}_0 \mathcal{S}^2$.
  and that the operator $\mathcal{S}$ commutes with $\mathcal{H}$ and $\mathcal{G}$.

How is this obvious? I have tried recreating these properties in R, but I cannot convince myself that $\mathcal{S} \mathcal{D}_0 = \mathcal{D}_0 \mathcal{S}^2$, or that S and the filters commute.  Can anyone help?

I am suspicious due to the claim that $\mathcal{R}_0$ is the inverse of $\mathcal{D}_0$. Decimation is not invertible. Once samples/entries are deleted, those values are forgotten. They cannot be reconstructed from what remains.

Let's say the input finite sequences are in $\mathbb{R}^{2N}$, numbered as $x = (x[0],x[1],\ldots,x[2N-1])$. Then $\mathcal{D}_0:\mathbb{R}^{2N}\to\mathbb{R}^N$, and $\mathcal{S}$ is overloaded, so that it maps from $\mathbb{R}^{2N}$ to $\mathbb{R}^{2N}$ and from $\mathbb{R}^{N}$ to $\mathbb{R}^{N}$:
\begin{eqnarray}
(\mathcal{S}x)[n] &=& x[n+1~\textrm{mod}~2N]~~\textrm{for}~~x\in\mathbb{R}^{2N},\\
(\mathcal{S}y)[n] &=& y[n+1~\textrm{mod}~N]~~\textrm{for}~~y\in\mathbb{R}^{N}.
\end{eqnarray}

Let $u = \mathcal{D}_0x\in\mathbb{R}^{N}$:

\begin{equation}
u[n] ~=~ (\mathcal{D}_0x)[n] ~=~ x[2n],
\end{equation}
where $n$ runs from 0 to $N-1$.

Then $\mathcal{S}u = \mathcal{S}\mathcal{D}_0x\in\mathbb{R}^{N}$:

\begin{equation}
\begin{split}
(\mathcal{S}u)[n] &=~ u[n+1 ~\textrm{mod}~N]\\
 &=~ x[2\times((n+1)~\textrm{mod}~N)]\\
 &=~ x[2n+2~\textrm{mod}~2N]
\end{split}
\end{equation}

Now let $v = \mathcal{S}^2x\in\mathbb{R}^{2N}$:

\begin{equation}
v[n] ~=~ (\mathcal{S}^2x)[n] ~=~ x[n+2~\textrm{mod}~2N]
\end{equation}

Then $\mathcal{D}_0v = \mathcal{D}_0\mathcal{S}^2x\in\mathbb{R}^{N}$:

\begin{equation}
\begin{split}
(\mathcal{D}_0v)[n] &=~ v[2n]\\
 &=~ x[(2n)+2~\textrm{mod}~2N]\\
&=~ x[2n+2~\textrm{mod}~2N]
\end{split}
\end{equation}

Hence
\begin{equation}
(\mathcal{D}_0\mathcal{S}^2x)[n] ~=~ x[2n+2~\textrm{mod}~2N] ~=~ (\mathcal{S}\mathcal{D}_0x)[n]
\end{equation}
for $0\leq n < N$ (because all results are in $\mathbb{R}^{N}$), so $\mathcal{D}_0\mathcal{S}^2$ and $\mathcal{S}\mathcal{D}_0$ are the same.


Convolution of delta function in discrete time with parameter If I have a delta function of the following type

$$ n\delta(5n-20) * 3^{5n}u[n] $$ How can I calculate this convolution?

I have thought of the property $$ \delta[n-n_0] *x[n] = x[n-n_0] $$

Then if I say $$ k=5n $$ and $$ t[n]= 3^{5n}u[n]  $$

It should be that the convolution is equal to :

$$ nt[\frac{n}{5}-20] $$ and replacing t $$ n 3^{5(\frac{n}{5}-20)}u[\frac{n}{5}-20]$$

Is this right? Or am I wrong somewhere?

You can use the following argumentation to find the result. The discrete time unit-sample function $\delta[n]$ has the following property for integer $M$:
$$ \delta[Mn] = \delta[n] $$ and more generally you can conlcude that for integer $M$ and $d$ we have
$$ \delta[M(n-d)] = \delta[n-d] $$

Therefore you can replace $\delta[5n-20] = \delta[5(n-4)]$ with $\delta[n-4]$ and proceed as usual to find the result of the convolution as:

$$y[n] = n\delta(5n-20) \star 3^{5n}u[n] = n\delta(n-4) \star 3^{5n}u[n] = 4 \cdot 3^{5(n-4)}u[n-4] $$

where in this last line we have used the sifting and shifting properties of the impulse function.


Deconvolution of signal which is broadened by laser spot size I have a signal that I measured by optical means using a focused laser. The laser scans over a microscopic feature (in step size of 0.33 micron) and I have plotted the result of the measurement as a function of distance (x scale is microns):

Signal image:

I am interested in the width of the feature I am measuring, however in this case the width seems on the same order (or even much smaller) than the laser spot size. This results in the signal being a Gaussian shape which could potentially just be the shape of the laser spot. But I cannot know that for sure, because I'm trying to measure the width of this feature. It is entirely plausible that the feature itself gives a Gaussian response, I expect it is not an extremely sharp edge. But I'm trying to determine how sharp the edge is exactly, and the laser spot size is now limiting.

Assuming the feature is extremely thin and sharp, then the measurement would indeed yield this Gaussian shape. However the feature is definitely not infinitely sharp and I believe I can measure more accurately what its shape is by using deconvolution.

I have tried to measure the laser spot size (this would be the system response function, I guess?) by scanning over a feature that IS definitely 100% sharp. If everything was perfect the signal from that would be a step function, but of course due to the large laser spot size I actually get a sort of error function shape from which I believe I can determine the laser spot size:

Sharp edge signal image:

Unless I'm mistaken the parameter C in my fit gives me the 1/e^2 beam width (0.785 micron). I think I can use this to form the system response function, assuming that my laser beam is Gaussian (which is probably not true but that's the best I can do).

The original signal I measured should be a convolution of the actual feature and this Gaussian system response, so I should also be able to deconvolve the signal with this same Gaussian system response and therefore obtain a better measurement of the feature, right?

How would I go about doing this? I have the original signal data (image 1, distance vs signal) and I can assume a Gaussian spot for my laser as the system response function with the known size obtained from the measurement of the perfectly sharp edge (image 2), as far as I understand (not much more than this...) this should be enough to do a deconvolution?

There was a similar question here a while back with a conceptually similar problem.
Link to the earlier question.

That one was about scanning a target with an electron beam.
The question there was how to determine the true shape of the beam if you use a known target.  You would seem to have a beam with a known shape.
You've got it a bit easier, since you are dealing with a single dimension rather than a 2D shape.

The easiest way to start would be to look at examples for the Matlab deconv method.

Basically, you need two vectors - one that represents the beam, and one that represents the object as scanned with the beam.  Feed them into the deconv method, and you get a vector that represents just the object.

One way to go at it would be to create a vector for an object of known size - say you have a sample of exactly 5 microns width.  You make a vector representing say 10 microns, and you put a 5 micron wide step in the middle of it.  Now you scan for real a 5 micron wide object, and deconvolve that vector with your generated 5 micron vector.  Now, you have the beam vector to use when measuring unknown objects.

If you don't have Matlab, you could use Gnu Octave which works almost just like Matlab - many Matlab examples Just Work in Octave.

Example Octave code:

calibration_target=[1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1];
simulated_beam=[1,1,2,2,2,1,1];
calibration_scan=conv(simulated_beam,calibration_target,"full");
[recovered_beam_shape,beam_remainder]=deconv(calibration_scan,calibration_target);

test_item=[0,0,0,1,1,0,0,0];
simulated_testscan=conv(simulated_beam,test_item,"full");
corrected_test_item_scan=deconv(simulated_testscan,recovered_beam_shape);

subplot(4,2,1);
plot(calibration_target);
title("calibration target");

subplot(4,2,2);
plot(calibration_scan);
title("simulated scan of target with simulated beam");

subplot(4,2,3);
plot(simulated_beam);
title("simulated beam");

subplot(4,2,4);
plot(recovered_beam_shape);
title("recovered beam shape");

subplot(4,2,5);
plot(simulated_testscan);
title("simulated testscan");

subplot(4,2,6);
plot(corrected_test_item_scan);
title("corrected test item scan");

subplot(4,2,7);
plot(test_item);
title(" test item true shape");

Result:

Added to answer questions from the comments.

Convolution and deconvolution work on uniformly spaced vectors.  All you need are the intensity values - drop your x coordinates.  If your data points aren't uniformly spaced, then you've got a problem.  If you interpolate your points to get uniform spacing, then you will probably introduce some fuzziness.  I don't know how much of a bad effect that would have on the deconvolution.  As for determining the vector to correct for the beam, I would scan a target with a known shape and sharp edges.  Then I would deconvolve the scanned data with an ideal image of the target - basically a vector as wide as the scanned object, with values of zero for the areas outside of the target, and values of 1 inside the target.  For example, a target 5 microns wide scanned with 1 micron stept would be [0,0,0,0,0,1,1,1,1,1,0,0,0,0,0]  Since you did a scan with a known object in order to generate your gaussian, I assume this is possible.
The "size" shouldn't matter - assuming you mean the amplitude or intensity.  I would scale the generated data to make 1 represent the scanned area at maximum intensity.  If you meant the number of elements, then that would matter.  You need to overscan by enough to be sure that your beam actually leaves the target to each side.  If you don't overscan enough, then you won't be able to get a true image of your target.


Why is cos(n/6) aperiodic? This is a very common example in most Signal Processing books I have come across.

x(n) = cos($\frac{n}{6}$) is a non-periodic discrete signal because it doesn't satisfy the periodicity condition for discrete time signals i.e, it is not of the form 2$\pi$($\frac{m}{N}$).

My question is :

the coefficient of n, i.e, $\Omega_0$=$\frac{1}{6}$ here can also be expressed as $\frac{1}{6}$ = $\frac{1}{6}$ * $\frac{2\pi}{2\pi}$ = 2$\pi$$\frac{1}{12\pi}$

Now, substituting for $\pi$ = $\frac{22}{7}$ in above, we get 2$\pi$$\frac{7}{12*22}$. So, $\frac{1}{6}$ can be written as 2$\pi$($\frac{7}{264}$), which is in the form 2$\pi$($\frac{m}{N}$) with a period N=264.

I'm sure I'm missing something which may be obvious but it would be of great help if someone could point it out and explain.

The problem with your reasoning is that $\pi \ne \frac{22}{7}$; $\pi$ is an irrational number. There is no period $N$ for which $x[n] = x[n+N] \ \forall \ n \in \mathbb{Z}$. Hence, the sequence is not periodic.


difference between Gabor and log-Gabor function I am reading a paper using log-gabor filters for feature detection. I was thinking about the difference between Gabor filters and log-gabor filters. Can anyone tell me the difference(s), and a way to implement them? Thanks in advance.

The 1D gabor filter has the following form in the frequency domain:

$$G_{b(\sigma,\omega_0)}(\omega) = \text{exp}\left(-\frac{\sigma^2}{2}(\omega - \omega_0)^2\right)$$

The 1D log-gabor filter is:

$$G_{l(\sigma,\omega_0)}(\omega) = \text{exp}\left(-\frac{\ln^2(\omega/\omega_0)}{2\ln^2(\sigma)}\right)$$

Log-gabor filters are used because they have 0 DC component for arbitrary large bandwidth, and size distribution of features in an image is often logarithmic.

The excellent paper On the Choice of Band-Pass Quadrature Filters by Djamal Boukerroui, J. Alison Noble and Michael Brady explains further differences for these two filters as well as others.

Here is some MATLAB code for the 2D version:

[rows, cols] = size(I);

% FFT mesh
[ux, uy] = meshgrid(([1:cols]-(fix(cols/2)+1))/(cols-mod(cols,2)), ...
    ([1:rows]-(fix(rows/2)+1))/(rows-mod(rows,2)));
ux = ifftshift(ux);   % Quadrant shift to put 0 frequency at the corners
uy = ifftshift(uy);

% Convert to polar coordinates
th = atan2(uy,ux);
r = sqrt(ux.^2 + uy.^2);

% Create spectrum
filterFFT = 1.0/wavelength;
filterFFT = exp((-(log(r/filterFFT)).^2) / (2 * log(sigma)^2));
filterFFT(1,1) = 0;

% Filter image
I_filtered = real(ifft2(fft2(I) .* filterFFT));

Advantage of 0 DC

The advantage of having 0 DC component is that the response of the filter doesn't depend on the mean value of the signal. It also means you don't have an infinite impulse response. The Gabor filter has a small non-zero component around DC that is usually removed.

For feature detection, say you designed a filter kernel to match a step (1D) or edge (2D). If the filter had a DC component then the response would vary with the mean image level. e.g. consider the sobel edge detector, and the same kernel with increased values.

I = double(imread('the_image.jpg'));
I = I(:,:,1);
f1 = [[-1 0 1];...
      [-2 0 2];...
      [-1 0 1]];
f2 = [[0 3 4];...
      [1 3 5];...
      [0 3 4]];
I1 = conv2(I,f1,'same'); images(I1); colorbar; pause;
I2 = conv2(I,f2,'same'); images(I2); colorbar; pause;
I3 = conv2(I+100,f1,'same'); images(I3); colorbar; pause;
I4 = conv2(I+100,f2,'same'); images(I4); colorbar; pause;

In this code f2 = f1 + 3. Only convolving with f1, which has 0 DC, gives the same result when the mean value of I is changed. The image responses 'look' the same, but only the 0 DC one is useful for feature detection.


Why is the concept of a "state covariance matrix" necessary in estimation? I'm currently taking a course in optimal estimation (and it's still very early in the course).

Much of our work is based around the idea of a measurement model

$y=Hx + v$

This model assumes our measurement $y$ is a function of some "state" and measurement noise.

Now when attempting to estimate the state $x$ in a stochastic least squares sense, one of the goals is to minimize the trace of the state covariance matrix $P$.

My question revolves around the idea of a state covariance matrix. Why exactly does it make sense to assume such a thing exists? By doing so we making the assertion that our "state" is a random variable.

Let's say we are trying to estimate the resistance of an unknown resistor by inserting a known voltage and measuring the current. In this case, our "state" is the resistance. Why would it make sense to assume our resistance is a random variable?

Do we assume the state is a random variable since there might be uncertainty in our model?

A resistance isn’t a particularly dynamic state, should be an unknown constant, but you might have a bin of resistors where they might vary. Taking one “randomly” out the bin makes it a random variable.  The bin of resistors will have a mean value, so perhaps that mean constitutes a state variable.  Perhaps someone starts putting resistors in the bin from a different source where they tend to run high, so the mean changes.

You don’t “have” to assume that the state is random if the model isn’t really random. Later you will find out in your class that states do not have to be stochastic. The KF equations will work for deterministic states. The measurement does need to have noise. If it doesn’t have measurement noise, you use a state observer, not a KF.

A random state is actually very general and modeling errors are often successfully considered noise.

Typically we use states because they encapsulate all the previous states. You use the current state to move to the next state. you don’t need to consider any past states to go to the next state. Some people will say a state has memory. If the state is stochastic, it has the Markov property. Your measurement depends only on the current state.

The idea of randomness has many interpretations. It is part of what makes the topic interesting.


Difference in signal Spectrum [
Here are two figures from the same signal, the 2nd one(B) is after digital predistortion at transmitter end. How can you describe it? What are the improvements and advantages of this technique?
Much appreciate your answers.

It is a frequency dependent gain to compensate attenuation of the high frequencies in the signal band (< 0.6Hz). As a side effect it also increases the background noise (>0.6 Hz)


Taking the FFT of a sinusoidal signal and going back I am a computer science student and want to do some stuff with audio data. I want to use the DFT to analyze and synthesize some sounds. Before going to the more complex stuff I experimented with basic tones and found some issues I do not understand - therefore I want to ask for help in this physics forum.

I created a single sinus wave function with an audioprogram of $440\textrm{ Hz}$. From this I take the DFT and want to recreate the original signal from this. While doing so, first I emit the negative frequencies, thus take only the first $n/2$ frequencies. $N$ is my blocksize (I divide the signal into blocks of length $n$). The sampling rate of my audio signal is $44100\textrm{ Hz}$.
So when I do this, and take the IFT from the resulting data, I can perfectly reconstruct my signal.

Now to my issues and questions:

When I set the phase to $0$ in the Fourier space (the imaginary part), with different $n$, I recover a sinus signal with some periodic noise.

I assume this is the phase difference which can be heard between the blocks?

When I set $n$ to $44100$ (the sampling frequency), I get complete noise.

Why is that?

Now I want to take only the strongest frequency amplitude wise (which in my opinion should perfectly work) - thus I set all other frequencies in the Fourier space to $0$ and then do the IFT. For different $n$ this kind of works, I get a sinus signal with some periodic noise.

Why the noise?
Moreover, the frequency of the resulting tone changes with $n$. Why is that?
With $n = 44100$ I get complete silence. Why?

When I set the phases to $0$ again I get the same results.

I hope it's somewhat interesting for you too. Could you explain to me these "phenomena"?

Part of the problem is that you are using sine waves.  The imaginary components of an FFT result contains all the information about pure sine waves (e.g. that are equivalent to a sin() function that starts with phase of zero at the start of the FFT aperture, and that are exactly integer periodic in aperture).  If you use only the real part of an FFT result, you can only see cosine waves (not sine waves, periodic in aperture), thus your silent result.

For shorter FFTs, the real component of the result might contain windowing artifacts if the frequency of your sine wave is between FFT result bins for the given FFT length (e.g. not integer periodic in aperture).  So, in synthesis, you are hearing an artifact, quantized in frequency.


Is the Laplace transform redundant? The Laplace transform is a generalization of the Fourier transform since the Fourier transform is the Laplace transform for $s = j\omega$ (i.e. $s$ is a pure imaginary number = zero real part of $s$).

  Reminder:

  Fourier transform: $X(\omega) = \int x(t) e^{-j\omega t} dt$

  Laplace transform: $X(s) = \int x(t) e^{-s t} dt$

Besides, a signal can be exactly reconstructed from its Fourier transform as well as its Laplace transform.

Since only a part of the Laplace transform is needed for the reconstruction (the part for which $\Re(s) = 0$), the rest of the Laplace transform ($\Re(s) \neq 0$) seems to be unuseful for the reconstruction...

Is it true?

Also, can the signal be reconstructed for another part of the Laplace transform (e.g. for $\Re(s)=5$ or $\Im(s)=9$)?

And what happens if we compute a Laplace transform of a signal, then changing only one point of the Laplace transform, and compute the inverse transform: do we come back to the original signal?

The Fourier and the Laplace transform obviously have many things in common. However, there are cases where only one of them can be used, or where it's more convenient to use one or the other.

First of all, even though in the definitions you simply replace $s$ by $j\omega$ or vice versa to go from one transform to the other, this cannot generally be done when given the Laplace transform $X_L(s)$ or the Fourier transform $X_F(j\omega)$ of a function. (I use different indices because the two functions can be different for the same time domain function). There are functions for which only the Laplace transform exists, e.g., $f(t)=e^{at}u(t)$, $a>0$, where $u(t)$ is the Heaviside step function. The reason is that the integral in the definition of the Laplace transform only converges for $\Re\{s\}>a$, which implies that the corresponding integral in the definition of the Fourier transform does not converge, i.e. the Fourier transform doesn't exist in this case.

There are functions for which both transforms exist, but $X_F(j\omega)\neq X_L(j\omega)$. One example is the function $f(t)=\sin(\omega_0t)u(t)$, for which the Fourier transform contains Dirac delta impulses.

Finally, there are also functions for which only the Fourier transform exists, but not the Laplace transform. This means that the integral in the definition of the Laplace transform only converges (in a specific sense) for $s=j\omega$, but for no other values of $s$. The Laplace transform is only said to exist if the integral converges in a half-plane or in a vertical strip of finite size of the complex $s$-plane. Such functions for which only the Fourier transform exists include complex exponentials and sinusoids ($-\infty<t<\infty$), and impulse responses of ideal brick-wall filters, which are related to the sinc function. So, e.g., the functions $f(t)=\sin(\omega_0 t)$ or $f(t)=\sin(\omega_ct)/\pi t$ do not have a Laplace transform but they do have a Fourier transform.

The Laplace transform can be a convenient tool for analyzing the behavior of linear time-invariant (LTI) systems by considering their transfer function, which is the Laplace transform of their impulse response. The poles and zeros of the transfer function in the complex $s$-plane conveniently characterize many system properties and are useful for an intuitive understanding of the system's behavior. Furthermore, the unilateral Laplace transform is very useful for analyzing LTI systems with non-zero initial conditions. The Fourier transform is a useful tool for analyzing ideal (non-causal, unstable) systems, such as ideal low pass or band pass filters.

Also have a look at this answer to a related question.


Are there any pictures of a sample accurate spectrogram on the web or in books? I programmed a sample accurate spectrogram some time ago and i still have never seen a picture of another one. It's interesting maths, i was told it's trivial, so why has no one else posted research on sample accurate sgrams on the internet? Can you find one. Here is a pic of my version:

this is the same program without the zoom and in 2d Mode:

Producing accurate time-frequency analysis is an area of research. Some pointers are : Wigner-Ville transformation, re-allocation methods, wavelets.


Tensor in signal processing I saw a lot of use of tensor in signal processing. What is the intuition behind it? Is it simply a common representation for audio (1D), image (2D) and video (3D) signals?

This has little to do with intuition. Tensors are rigorously defined mathematical objects, and in general simple arrays don't qualify as tensors. Specifically, signals are not tensors.

In the language of mathematics, and there the field of differential geometry, a tensor is a linear function with several arguments. It is therefore also called a multi-linear form. It is used to describe properties of manifolds at a single point.

A point on a (smooth) manifold comes with a tangent space that contains all possible "direction" at that point, where a direction is really just a tangent vector at that point. This tangent space at a single point is a vector space and it has the same dimension as the manifold.

The tangent space has a natural so called dual space, which contains the linear functions that map tangent space vectors to real numbers. This is the cotangent space, and it's also a vector space of the dimension of the manifold.

A tensor is now (roughly) a linear function that maps any number of dimensions between the tangent space, the cotangent space and the real numbers.

Tensors typically don't come as singles but in so called bundles, which are the disjoint union of the tensors of all points of the manifold. For example if you have a vector at every point of the manifold then the whole structure is called a vector bundle. Similarly, tensor bundles associate a tensor with every point on the manifold, usually in a smooth way.

In signal processing you can encounter tensors, but usually they are called differently. For example in video processing the image is your manifold (or rather, a function on a 2-dimensional manifold giving the brightness or color for each point) and the velocity field that describes the local motion is a rank 1 tensor field (or bundle) on that manifold.

Tensors also show up a lot in volume data processing, where tensors can describe local properties of the volume data. For example if you have a doppler ultrasonic tomography image, then the velocity data in each voxel is a tensor field as would be mechanical stress in tomographic material analysis.

For a 1-dimensional signals tensors often come in the form of derivative operators. For example if you have a signal $s(t)$, then the operator $\frac{\partial}{\partial t}$ is a tensor field of rank 1, or just a vector field. More generally, if you multiply that operator with a smooth function $g(t)$, then all possible vector fields on the manifold the signal lives on are of the form $g(t)\frac{\partial}{\partial t}$

To sum up, tensors are a description of differential properties of manifolds. If you want to understand them, you'll have to study differential geometry.


Causality of z-transform $a^nu[n+1]$ To preface, this is not a homework related question but purely for self-study purposes.

I'm try to do the analyse of z-transform of $a^nu[n+1]$. It is clearly a non-causal signal, I try to explain it by using the definition.

Based on the causality of z-transform definition: If x(n) is a infinite duration causal sequence, ROC is exterior of the circle with radius a. i.e. $|z| > |a|$.

I am not sure if my explanation is correct:
$$\mathcal{Z}\big\{a^nu[n+1]\big\} = a^{-1}z+\sum_{n=0}^{\infty} (a/z)^{n}$$

The ROC of this function in z-plane is $|z| > |a|$ except at $z = ∞$.

It does not obey the causality definition because the boundary exists ( $|a|<|z|< ∞ $), i.e. the ROC pattern is a ring, it is a non-causal signal.

Am I correct? I am kinda confused with the causality definition.

You are correct. If the region of convergence of a right-sided signal (like the one you have) does not include infinity, then the signal is not causal.


Generating band-limited NRZ signal What is the bandwidth of a NRZ signal (of a given frequency) and how would you directly synthesize a discrete version a band-limited NRZ signal at some higher sample rate?

The first concept that I think will be helpful to see is pulse shapes that repeat in time at consistent rates. If you repeat a fundamental pulse shape in time, the result in frequency will have the envelope of the Fourier Transform for that pulse but will only contain frequency components at multiples of the repetition rate (harmonics).  Thus the Fourier Transform becomes sampled in frequency at multiples of the repetition rate in time. For example see the cases below for NRZ signals with different duty cycles, knowing that the Fourier Transform of a square pulse is a Sinc function in frequency with the first null at $1/T$ where T is the width of the pulse.

Now consider what occurs if the NRZ data has a pulse width of 0.05 seconds like the first example but is a random pattern 101101110001... that doesn't repeat for 10 seconds - the envelope will still be the same as the first example with the first null at 20 Hz, but the harmonics will appear every 0.1 Hz, the spectrum will appear continuous. Extend this thought to random data that never repeats and the spectrum will indeed be continuous with the same shape.

So to answer your question the bandwidth of an NRZ pulse is infinite, but the bandwidth of any time limited pulse will be infinite if you review the Fourier Transform equations. However for practical considerations we can define bandwidth to be % of power within a certain frequency band, or ignoring spectral content that gets below a certain insignificant value.

With that in mind, the bandwidth of a Sinc pulse in frequency (NRZ signal in time) is still quite large, since the peaks of the Sinc only roll-off at rate f in frequency (or 20 dB/decade).  By shaping the pulse edges in time (slowing the transition rate), we can significantly limit the bandwidth.

A common approach is to reverse the case of the NRZ signal above, and use time-limited versions of the Sinc function in time as the fundamental pulse, in such a way that the null crossings of the Sinc are aligned with the centers of a subsequent pulses (for zero inter-symbol inteference while getting the frequency limiting abilities of the Sinc pulse). The longer you can make the Sinc pulse in duration, the tighter the bandwidth as in the extreme of an infinitely long Sinc in time would result in a brick-wall spectrum in frequency.

This is demonstrated clearly in this graphic below, for a muli-level PAM waveform (Or what it actually is, is the real part of a 16-QAM waveform). Regardless the fundamental pulse shape for the blue line is a rectangular NRZ pulse in time, and the blue spectrum shown is what I was explaining above for the case or random non-repeating data.

The red line is the "Raised-Cosine" pulse shape which is a Sinc pulse that has an additional taper given by a "raised cosine" (see additional graphic on that). But notably see how significantly more constrained the resulting red spectrum becomes and hence the pulse-shaped waveform is more spectrally efficient. We go through the extra effort to shape the pulse in order to conserve spectrum.


Perform frequency analysis on grouped pulses I have a system which consists of individual pulses grouped in trains.

The trains have a frequency of 10 Hz, with a timing precision of sub ns.
The pulses have a frequency of 2.2 MHz with a timing precision of sub ps. The train consists of 100-1000 pulses.

I would like to do a spectrogram in the range 0.1 - 100 Hz. Up to  5 Hz its trivial. I simply take the average of the train and still have the Nyquist theorem fullfilled. However I see a lot of signal which could come from aliasing. Is there a possibility to use the 5 order of magnitude faster signal from the individual pulses to gain anything useful here? I have a maximum of 25d of data. But I would ideally like to do a spectrogram, meaning it would be great to stay within about an hour of data (e.g 36k trains)

The first step is to low pass filter the signal with 100Hz cutoff and then resample everything at 200Hz. This doesn't change the spectrum in the 0-100Hz range. Therefore, it's not possible to gain anything of use from the faster signals.

The 0-100Hz spectrum will be the same as if there were a single 45-450μs pulse occuring at 10Hz. There will peaks at multiples of 10Hz (harmonics) and the maximum will be the 10Hz peak itself.


Low MIPS Video Encoder I am looking for a low MIPS, low compression video encoder. This is for a 10fps, VGA type quality compression. What are my options? I need to be able to do this compression using a 150MHz ARM M4 CPU with Floating point support. (STM32F4)..

My idea is to push this compressed data out of the CPU with a parallel bus, there will be no processing done on the data. Compression rate wise, as much as possible, I just want to see the boundaries. This is for a low cost CCTV application, I like to see what I can achieve with a 5USD cpu and lots of transmission band with vs. a 30USD encoder with low data trasnmission bw.

10fps, VGA will generate roughly 25Mbit/sec data. This is fairly high data rate for anything out there. If I can bring this down to 5Mbit/sec, I think I can build a very low cost CCTV system. Once I get the data to the base, I can re-encode the data, therefore it doesn't matter what the compression mechanism is, as long as it is not very lossy.

Monochromatic video is what is necessary at moment over color.

Update

This CPU has 120MHz allocated to this task.
Memory interface is 16 bit, therefore external memory write/read is slower compared to internal memory.
Internal memory is 120KByte and it has 32 bit access fast access. Both cases the memory is accessed via AHB bus which we should assume 60MHz as the clock frequency.
Following data flow is expected:

Camera ->DMA ->External Memory (no cpu involvement)
External Memory -> CPU ->Compression ->Internal Memory
Internal Memory -> DMA -> Data Bus -> External Device

CPU will only read a chunk of data compress and write to its internal memory (the compressed data), later start a DMA transfer.

Your first problem is memory - even the largest version of that processor doesn't have that much internal SRAM (in video terms - it's a stonking embedded processor compared to many though!) - a VGA frame is 300kB, vs the 100ishkB in the processor.  That means you'll have to process it as it comes in - say in 8 or 16 line chunks, which makes it much more of a "hard-realtime" problem as the deadlines have less slack in them.

It's not clear how you'll be capturing data, but I assume that you'll have some sort of DMA from either the ADC or some digital port, otherwise you'll be spending half your time shuffling data around!

Aiming for lowish compression ratios, you might just want to look at encoding the deltas between pixels along a line or a 3x3 square and then arithmetic code them.

See also Simple, streaming, lossless image compression - although I was asking for explicitly lossless compression, there may be some ideas for you there.

As a datapoint, one of my colleagues implemented a 320x240 8-bit mono JPG compression in a 60MHz 32-bit micro with floating point unit and a small cache.  He used some reference code (ie optimised for readability, not performance) and got 5 fps fairly easily.  Compression ratios were of the order of 10x IIRC. The image capture was done by external hardware bus-mastering the data into the micro's external RAM.


Voice Recognition - MFCC I have been researching about Speech Recognition and I have decided to take the MFCC approach to solving this problem of detecting whether someone is saying either "Yes" or "No". So as mentioned before (My steps so far):

Read in Audio File
Split the Audio Signal into blocks (600 samples, 30msec long)
Strip the blocks that do not warrant consideration (Total Energy / Zero-crossing)

So I am going to construct the MFCC based on this paper and it has the following steps:

Pre–emphasis
Framing
Hamming Windowing
FFT
Mel Filter Bank Processing
Discrete Cosine Transform
Delta Energy and Delta Spectrum

This makes sense to me (Kind of) and I am going to research into each of these steps. BUT should I perform the MFCC on the resulting blocks that I have already done with (steps 1, 2, 3) at the top of this question, or, should I not carry out these steps and just start from the beginning and compute the MFCC and will I still be able to implement a Hidden Markov Model?

The other question is, if I split the signal into "Frames" (2D vector) will the resulting MFCC be a 2D vector, or a 1D vector?

Hope someone can help :)!

Regarding your existing steps:

Do you have some overlap between adjacent blocks? It is common in feature extraction systems to have some overlap between adjacent blocks, so that a short transient event can be correctly captured (if it is right at the end of a block, it will be right in the middle of the next one).
Your idea of discarding blocks is dangerous. You can do this to roughly identify the beginning and end of a speech segment; but silence recognition can also be built into the recognition model.

Regarding the data size:

Your signal is split into frames. For each frame, you compute a 1-D MFCC vector. So in the end, you have a sequence of 1-D MFCC vectors. This is the data on which you train your HMMs.


How to minimize peak-to-peak amplitude of periodic waveform through FFT phase shifts? For one of my music projects, I'm playing back periodic audio signals (by looping single periods). Unfortunately, one of my waveforms sounds too quiet (even at maximum volume).

I'm trying to use FFT to obtain harmonic strengths, then phase-shift each harmonic to minimize the peak-to-peak amplitude, and then amplify the signal.

Is there any algorithm to find the optimal phase shifts, or should I use scipy.optimize.basinhopping (or generating thousands of random phases) for an approximate result?

There are algorithms for determining low-peak-factor signals, for example:

  "Synthesis of low-peak-factor signals and binary sequences with low
  autocorrelation (Corresp.)", M. Schroeder, 1970

I have used this for flat-spectra signals to minimise peak to peak amplitude.

However I think this may be an XY problem. Do you have a limit on how much distortion you can introduce to the signals? Have you tried using a dynamic range compression algorithm? There are limiting algorithms which can increase the RMS amplitude significantly without introducing significant audio artifacts.


Sampling and recovery using matlab I'm trying to write a program using MATLAB that samples and recovers a given signal with different sampling frequencies

The given signal is $x(t) = \sin (2\pi 200 t)$ and sampling frequencies are $200$,$300$,$400$ and $500$ Hz.

I used sinc interpolation formula.

It is not working,and I would like to know why.

Can you help me?

The code is:

 f=200;
 fs=[200,300,400,500]; % sampling frequencies
for i = 1:length(fs)

 Ts = 1/fs(i); % Sampling Period
 t = 0:1e-4:6/f; % time vector for continous plot
 w = f*2*pi;
 xa = sin(w*t); % signal in continuous form

 t2=(0:Ts:6/f);  % time vector for sampling
 xb = sin(w*t2); % sampled signal

 figure
 subplot(311);
 plot(t,xa) % continous signal plot
 str = ['Signal x(t) = (' num2str(f) ' Hz)', ', Fs = ' num2str(fs(i)) ' Hz'];
 title(str);
 ylabel('Amplitude');
 xlabel('Time(s)');

 subplot(312);

 stem(t2,xb); % sampled signal plot
 title('Sampled signal');
 ylabel('Amplitude')
 xlabel('Samples');

 xr = zeros(length(t));
 for k=1:length(t)
 for p=1:length(xb)
 xr(k)=xr(k) + xb(p)*sinc(fs(i)*(k-t2(p)));
 end
 end
 subplot(313);
 plot(t,xr);
 title('Recovered signal');
 ylabel('Amplitude');
 xlabel('Time(s)');
end

This line is wrong:

xr(k)=xr(k) + xb(p)*sinc(fs(i)*(k-t2(p)));

The correct one is:

xr(k)=xr(k) + xb(p)*sinc( ((k-1)*Tr-(p-1)*Ts)/Ts );

Where $Ts$ is the original sampling period and $Tr$ is the reconstruction sampling period. In your code $Ts$ is given by $T_s = 1/500$ for sampling frequency of $500$ Hz. And $T_r = 10^{-4}$ set by your continuous plot period.

Note that you will still have slight errors due to the fact that the input signal is not being perfectly bandlimited.

Here is the resulting plot for $Fs=500$ Hz, $f_0 = 200$ Hz and $T_r=10^{-4}$ seconds:

it's ok now.


Reconstructing/interpolating small regions of a bandlimited signal by taking the fewest possible samples I have a signal which is bandlimited and can be sampled at arbitrary continuous positions. The value at any position is given by an expensive computation. I need to do some further computation on various arbitrary small regions of this signal.

My plan is to take advantage of the known bandlimiting to sample the signal at or slightly above the Nyquist rate within these regions that I need to process, and then do my computation on a reconstructed approximation of the region. This reconstruction needs to be in a form that allows for evaluation at arbitrary time t within the reconstructed range. Because of the expense of computing each sample, I need to minimise the number of samples required to reconstruct each region. The regions are sized such that they are between 8 and 16 samples wide at the Nyquist rate.

The thing that makes this different from standard usages of sinc interpolation is the small number of samples. That's the focus of this question. It's typical for sinc based reconstructions to experience problems near the edges of the signals, but when a signal is only 8 samples wide, most of it is edge.

My question concerns how to optimally reconstruct these small regions with minimal samples. I particularly would like to know whether it's possible to achieve a reasonable reconstruction with all samples being within the region I'm reconstructing. I suspect that the answer is no, because there's insufficient information to reconstruct towards the edges of the region, but I would like to have this confirmed. My experiments with various windowed sinc methods support this suspicion.

Assuming that I can't rely on only samples within the region, is there anything better I can do than taking a few extra samples on either side of the region and then continuing to use basic sinc methods?

I don't have a firm definition of a good enough reconstruction, but for the purposes of the question, let's say that at all points along the reconstruction, the error should be no more than 10%.

I have some very short signals in the range of 8 to 16 samples. These represent a bandlimited signal, sampled at or slightly above the Nyquist rate.

Nope. A signal can't be limited in time and in frequency at the same time. If it's very short, than chances are the bandwidth is a lot higher than you think it is and that you've already picked up some significant aliasing in the sampling.

  I've experimented with sinc interpolation with a few different windowing methods, and also tried numerically optimising a kernel to achieve the interpolation, without success.

Can you elaborate what your problem is? If it's sampled correctly, sinc interpolation typically does a good job, although it does destroy any causality and will extent the signal in time. But that goes back to my first point: a band-limited signal can not be limited in time as well, so it can't be causal either.


One sided frequency spectrum (Matlab vs. Origin) There are a lot of queries on fft frequency all over the web. I guess the following point not discussed anywhere explicitly. Hope someone can provide an insight here.

If we have and even number of data points, N=10, the fft complex output arranges the data as

fft = [c0, c1, c2, c3, c4, c-5, c-4, c-3, c-2, c-1],

where the complex numbers corresponding to positive and negative frequencies. I read somewhere that MATLAB calculates the negative coefficient first, hence we have c-5 but not c5. The author did not explain the reason.

Point no. 1, that the values are not symmetric, there is an extra c-5.

When we wish to make two-sided frequency spectrum, the frequency axis ranges from [-(N/2): (N/2)-1]*Fs/ N. Fs is the sampling rate, N is the number of even data points.

If we wish to make a one-sided positive frequency spectrum, should we choose

A) [0:(N/2)]*Fs/N and ignore the fact the we are using the values corresponding to the negative frequency axis, given that the data is a real number and it is just a mirror image.

B) [0: (N/2)-1]*Fs/N represents the true positive frequency axis?

If Fs= 250 Hz, the true positive frequency axis will end at 124.9980 Hz

If we happen to choose the negative frequency axis values and ignore the frequency sign, the frequency axis ends at 125 Hz exactly.

The same data when plotted in OriginPro ends the frequency axis at 125 Hz when plotted single sidedly.

Which approach (A) or (B) is rigorously correct, and why would OriginPro plot the frequency upto 125 Hz?

The indexing is a matter of convention.  "Natural" is zero centered, "Computer implementation" is zero based.

We are having a big discussion right now over several questions on why your $c_{-5}$ should actually be considered as $(c_{-5} + $c_{5})/2$.  You'll find part of the discussion here Convergence of periodic sinc interpolation and the links to the rest are embedded within.

The "fluffy cloud" drawings show what happens when you consider the upper half (from a zero based perspective) as positive frequencies.  How to get Fourier coefficients to draw any shape using DFT?

If I have a circle, and I check your location every minute and it seems you have moved a quarter circle counter-clockwise each time, how do I know you haven't actually moved three quarters of the circle clockwise?

Answer:  You don't.  The convention decides which interpretation to use.


multiplications per second for FIR filter In my filter the fs=300MHz. and no. of coefficients is 31. I have following queries.
(1) what will be the number of multiplications per second?
(2) can multiplication be done in one clock cycle?
(3) why normally multiplication per second is calculated and not addition/substraction per second in the filter implementation?

To answer your questions:

The number of multiplies is the number of taps times the number of sample per second.  Given a sampling rate of 300MHz and 31 taps, you will have to do 9300 million multiplies per second.
As Paul R says, whether a multiply can be done in one cycle depends on the processor.  Some can, some can't, and some appear to by having a long pipeline - once the pipeline fills, the result of one multiply comes out per cycle.  You would have to find out what processor you are using to determine how many cycles it takes for each multiply.  Also note, this varies also with the type of number you are working with.  Integer multiplication is usually faster, but can cause you other difficulties (it is harder to do signal processing with integers.)
As Paul R says, really only the multiplies are of interest since most processors have a combined multiply and add instruction.  If you are working with a processor without that type of instruction, though, you will have to remember to count the adds when budgeting your clock cycles.  Since there are as many adds as multiplies when doing an FIR, knowing how many multiplies there are also tells you the number of adds.


Analytical way to generate an IIR that is the square root response of another IIR For simplicity take the transfer function of a second order IIR in $z$-domain:

$$
H(z) = \frac {b_0 + b_1 z^{-1} + b_2 z^{-2}} {a_0 + a_1 z^{-1} + a_2 z^{-2}}
$$

Is it possible to generate a new transfer function that yields the square root of the magnitude response of $H(z)$? Assuming the new transfer function is minimum phase. I know this would be possible using various numerical approaches such as global optimization or least squares, but has this specific problem been tackled or are their any approaches which would make the problem attainable without resorting to a heavy weight numerical approach? An approach that worked in the $s$-domain for an equivalent $H(s)$ would be equally useful.

The only time you'll be able to do this exactly is when you know that the poles and zeros of $H(z)$ are in complex conjugate pairs inside and outside of the unit circle.

Note that for systems with real-valued impulse responses, this means you need four of each: two inside and two outside the unit circle.

The keyword to search on is spectral factorization.

This means that if you can write $H(z)$ as:
$$
H(z) = G(z) G(z^{-1})
$$
where $G(z)$ is a minimum phase system, then $G$ just takes all the poles and zeros of $H$ that are inside the unit circle (and half those that are on the unit circle).

There exist approximate algorithms for finding a solution, but I am not well-versed in those.


How determine the delay in my signal practically I am trying to practically determine filter group delay so that I can try to reduce or adjust for it. I have tried all kinds of methods but still can't get the right result. I know that theoretically to get the group delay you should differentiate the phase response with respect to frequency. How can I use this theoretical result to determine what delay will actually occur to my whole signal envelope? I ask because the theory does not seem to apply in practice in the way I have understood it.

I have tried to align filtered and unfiltered signal to see an actual shift but I do not seem to see any delay.
I have tried adding a spike (Ok which is zero frequency) just so that I can compare the position of the spike in the original and the filtered signal, but the spike does not move in the filtered signal.

Additionally, some filters have both negative and positive theoretical group delay. In those filters should I take the average across all frequencies to determine the effective group delay?

I have added an example below in MATLAB. I used a linear phase FIR filter (a comb filter) which should have a delay of 20 samples according to MATLAB across all frequencies. But I have looked at the plots with the filtered and original signal on the same axis but could not get 20 samples at several frequencies.

%% Filter
b=[1,zeros(1,39),-1];%y(n)=x(n)-x(n-40)
a=1;

subplot(3,1,1)
grpdelay(b,1)

% Simulate
Fs=1000;
t=0:1/Fs:(5-1/Fs);
wi=blackman(length(t))';

spike=zeros(1,length(t));
spike(300)=0.02;%Feature
spike(150)=0.02;%feature

x1=sin(2*pi*49*t).*wi+spike;
x2=sin(2*pi*25*t).*wi;
x3=sin(2*pi*2*t).*wi;

x=[x1,x2,x3];

y=filter(cb,1,x);
subplot(3,1,2)
plot([x',y'])
title('x vs filtered x');
legend({'x','x-filtered'})

% show spike
subplot(3,1,3)
plot([x',y'])
title('x vs filtered x, zoomed to see spikes');
legend({'x','x-filtered'})
xlim([130,350])

My overall question is how can I practically measure the effective signal envelop delay due to filtering? I would like to match this measured delay with the theoretically calculated delay.

The Least Mean Square solution to find the "channel" or response of the filter is provided by the following MATLAB/Octave Code using  the input to the filter as tx and the output of the filter as rx. For more details on how this works, see this post: Compensating Loudspeaker frequency response in an audio signal:

function coeff = channel(tx,rx,ntaps)
    % Determines channel coefficients using the Wiener-Hopf equations (LMS Solution)
    % TX = Transmitted (channel input) waveform, row vector, length must be >> ntaps
    % RX = Received (ch output) waveform, row vector, length must be >> ntaps
    % NTAPS = Number of taps for channel coefficients
    % Dan Boschen 1/13/2020

    tx= tx(:)';   % force row vector
    rx= rx(:)';   % force row vector
    depth = min(length(rx),length(tx));
    A=convmtx(rx(1:depth).',ntaps);
    R=A'*A;       % autocorrelation matrix
    X=[tx(1:depth) zeros(1,ntaps-1)].';
    ro=A'*X;      % cross correlation vector
    coeff=(inv(R)*ro);   %solution
end

The case the OP uses of a comb filter is one of the most challenging as it depends on signal energy at each frequency for a solution (this is specifically why linear equalizers, which this function is doing if you swap rx and tx, do not perform well in frequency selective channels and result in noise enhancement at the null locations).  Below the frequency response of the test filter as determined with MATLAB or Octave showing the multiple frequency nulls associated with such a comb filter:

b=[1,zeros(1,39),-1];
freqz(b,1,2^14)     % 2^14 samples to show nulls

MATLAB or Octave script to demonstrate use of above function and determine the delay between output and input:

%% Filter with OP's example
b=[1,zeros(1,39),-1];     % numerator coefficients
a = 1;                    % denominator coefficients

%% Generate signal using OP's code
Fs=1000;
t=0:1/Fs:(5-1/Fs);
wi=blackman(length(t))';
rn=+rand(1,length(t))*.2;
x1=sin(2*pi*13*t).*wi +rn;
x2=sin(2*pi*25*t).*wi +rn;
x3=sin(2*pi*2*t).*wi +rn;
x=[x1,x2,x3];

% Filter
y=filter(b,a,x);

%% Test filter estimation
cf=channel(x,y,61);

%compare original and estimated channel
subplot(2,1,1)
stem(b)
title("Actual Channel Response")
xlabel("Sample Number")
subplot(2,1,2)

stem(cf)
title("Estimated Channel Response")
xlabel("Sample Number")

We could have used exactly 41 taps when calling the function channel and it would have resolved the solution properly; however it is best practice to start with a much longer filter length, evaluate the result and then reduce the taps accordingly. In actual practice under noise conditions using more taps than necessary will lead to noise enhancement, so it is good to do the final solution with the minimum taps necessary to capture the dominant tap weights.

Observe with the groupdelay plot using MATLAB and Octaves grpdelay command the issue with not being able to resolve the delay where no signal gets through the filter (it would be difficult for anything to determine the delay of a single tone at one of these frequencies that is nulled by the filter!) but it is able to accurately determined the delay where signal energy exists. Similarly the waveform itself must have energy at all frequencies where we are seeking a solution. The spectral density of the OP's test waveform was sufficiently spread out over all frequencies to be suitable for this purpose (this is the reason why pseudo-random waveform make for good "channel sounding" patterns).

This plot is for comparison to the OP's plot showing that the group delay for this filter is 20 samples.

Here is another test case showing how this works with a more reasonable channel (no deep nulls) with the following coefficients and frequency response below that:

b = [0.2 .4 -.3 .4 .3 .1];

The solution is indistinguishable between actual and estimate, so I added noise to make it more interesting, using same x and y from code above:

% add noise
noise = 0.351*randn(1,length(y));
yn = y + noise;
snr = 20*log10(std(yn)./std(noise));
%% Test filter estimation
cf=channel(x,yn,10);


Graph Fourier transform: the adjoint notation for the eigenbasis matrix I already asked this question here but there is no response. I'd like to ask this question in signal processing domain.

  It is well-known that for a real symmetric matrix $L$ (here, graph
  Laplacian) one can write the eigenvalue decomposition as

  $$ L = U \Lambda U^{\mathsf T}, $$ where $U$ is a real eigenvector
  matrix. Moreover, in graph signal processing papers, including the
  great paper by Shuman et al. (cf. page 4), the adjoint (complex
  conjugate) of $U$ is used to define the graph Fourier transform
  $\mathcal{F}_{G}$ as  $$ \hat{x} = \mathcal{F}_{G} x = U^{*}x, $$
  where $x$ is the signal in vector form and $U^{*}$ is the complex
  conjugate of $U$.

I am curious to know is there any specific reason for using the notation of complex conjugate since $U$ is real?

You terminology is something I am not used to.

$$ \hat{x} = \mathcal{F}_{G} x $$

If by $\hat{x}$ you mean the DFT (complex except very special cases), $x$ is your input (real or complex), and $\mathcal{F}_{G}$ is the matrix composed of the sinusoidal basis vectors, then your assertion that $U$ is real is dubious.  $\mathcal{F}_{G}$ is most certainly not.  With the proper normalizattion:

$$ \mathcal{F}_{G}^{-1} = \mathcal{F}_{G}^{*} $$

Now, for a known $x$ where $\hat{x}$ is known to be complex, how can $U$ be real?

Since they are using the complex conjugate, I would say that U is expected to also be capable of being complex.  I have not read the paper so I don't know where you picked up that assumption.


making pink noise (1/f) using list of frequencies making pink noise (1/f) using list of frequencies

I would like to see what type of noise I would get if I used just the frequency in my voice.  I created a matlab array using fft to get the
[frequency,amplitude,phase] to reproduce my vocal signal.

I would like to take this file/data and use it to create pink noise (1/f).  Of course when I use 1/f for the frequency the numbers become very small does anyone have any ideas how to use my own vocal frequencies that I get from doing a fft in matlab from my voice to create pink noise (1/f).

Thanks

Pink noise, as you know, has contributing amplitudes inversely proportional to each corresponding frequency. It appears more sinusoidal than white noise (constant power spectrum) because the low frequency components have large amplitudes. The high frequency components do not contribute that much to the overall noise amplitude.

To recreate pink noise in MATLAB, I would use the following technique:

1.) Generate a square wave (ones(n) will create a square wave of width n).

2.) Take the Fourier Transform of the signal (fft(x)). This will give you the sinc function. The DC component will be 1.

3.) In a for loop, add noise to the amplitude of each frequency of the FFT sinc function. The range of values of the noise will be dependent on 1/f.

4.) Take the inverse Fourier Transform of the noisy FFT. Then subtract 1 to yield a pink noise vector of length n centered around $0$.


Can Gabor filter be used for detecting dents in cars? I am doing some research on Gabor filters for detecting dents in cars. I know Gabor filters have had widespread use for pattern recognition, fingerprint recognition, etc.

I have an image.

Using some code from the MathWorks File Exchange site, I got the following output.

which somehow isn't the output that one would expect. This isn't a good result.

The my script is as follows:

I = imread('dent.jpg');
I = rgb2gray(I);
[G, gabout] = gaborfilter1(I, 2, 4, 16, pi/2);
figure
imshow(uint8(gabout));

EDIT: Applying a different code to the following image :

Output image after different orientations of gabor filter:

how do i isolate this DENT which is got being detected properly?

This is an extremely difficult problem. I was part of a team that worked on it for several years, and having developed and supported other such applications for a long time I can say that dent detection is a particularly tricky problem, and much harder than it looks at first.

Having an algorithm work under lab conditions or on known images is one thing; developing a system that is accurate and robust for "natural" images such as cars seen on a parking lot would likely require a team working for several years. In addition to the core problem of creating the algorithm, there are numerous other engineering difficulties.

The sample code you tested isn't a bad start. If you could find the edges around the dark right side of the dent you could compare the edge map of the car with the ding to an edge map of a known good car imaged from the same angle with the same lighting. Controlling the lighting will help quite a bit.

Problems to consider include the following:

Lighting (much more difficult than it would first seem)
Expected 3D surface of assembled outer panel (e.g. from CAD data)
Criteria characterizing a dent: area, depth, profile, etc.
Criteria for false negatives and false positives
Means to save dent data and/or map dents onto a model of the car (or butterfly layout)
Methodology and device to measure "true" dent characteristics: depth, area, etc.
Extensive database of dents from a random sampling of vehicles
Dealing with different paint colors and finishes

1. Lighting
As Martin B noted correctly above, correct lighting is critical for this problem. Even with good structured lighting, you're going to have great difficulty detecting small dents near feature lines, gaps between panels, handles, and so on.

The Wikipedia entry for structured lighting is a bit thin, but it's a start for understanding the principle:
http://en.wikipedia.org/wiki/Structured_light

Light stripes can be used to detect in-dings (dents) and out-dings (pimples). To see a ding, you'll need relative motion between the light source and the car. Either the light + camera move together relative to the car, or the car moves past the light + camera.

Although in-dings and out-dings have characteristic appearances when seen at the edge of a light stripe, the detectability of a given dent depends on the size and depth of the dent relative to the width of the light stripe. A car's curvature is complex, so it's quite difficult to present a consistent light stripe to a camera. As the light stripe moves across the car body, the curvature and even the intensity of the light stripe will vary.

One partial solution is to ensure that the camera and light stripe always at a consistent angle relative to the normal (the 3D perpendicular) of the portion of the surface being inspected. In practice a robot would be required to move the camera accurately relative to the body surface. Moving the robot accurately requires knowledge of the pose (position and 3D angles) of the car body, which is a nasty problem by itself.

For any inspection for automotive applications, you need to completely control the lighting. That means not only placing lights of your choice at known locations, but also blocking all other light. This will mean a fairly large enclosure. Since the car's panels are curved outward (almost like a spherical surface), they'll reflect light from sources all around them. To greatly simplify this problem, you could use a high frequency flourescent bar inside an enclosure shrouded with black velvet. Quite often it's necessary to go to extremes like that for inspection applications.

2. 3D surface
A car's outer surface is composed of complex curves. In order to know whether a suspicious spot is a ding, you have to compare that spot to known features of the car. That means you would need to match the 2D image from a camera to a 3D model viewed at a certain angle. This is not a problem solved quickly, and it's difficult enough to do well that some companies specialize in it.

3. Defect characterization
For academic research or lab testing it may be sufficient to develop an algorithm that shows promise or improves on an existing method. To properly solve this problem for real commercial or industrial use, you need to have a highly detailed specification for the size dents you'd like to detect.

When we tackled this problem, there were no reasonable industry or national standards for dents (3D deformations). That is, there was no agreed-upon technique to characterize a dent by its area, depth, and shape. We just had samples that industry experts agreed were bad, not too bad, and marginal in terms of severity. Defining the "depth" of a ding is tricky, too, since a ding is a 3D indentation in (typically) a 3D surface curving outward.

Bigger dings are easier to detect, but they're also less common. An experienced auto worker can scan a car body quickly--much more quickly than an untrained observer--and find shallow dings the size of your pinky finger quickly. To justify the cost of an automated system, you would likely have to match an experienced observer's ability.

4. Criteria for detection errors
Early on you should set criteria for acceptable false negatives and false positives. Even if you're just studying this problem as an R & D project and don't intend to develop a product, try to define your detection criteria.

false negative: dent present, but not detected
false positive: unblemished area identified as a dent

There's typically a tradeoff: increase sensitivity and you'll find more dings (decrease false negatives), but you'll also find more dings that aren't there (increase false positives). It is quite easy to convince oneself that an algorithm performs better than it actually does: our natural bias is to notice defects detected by the algorithm and explain away those it hasn't detected. Conduct blind, automated tests. If possible, have someone else measure the dings and assign severity so that you don't know what the true measurements are.

5. Save data and/or map it
A dent is characterized by its severity and its location on the car body. To know its location, you must solve the 2D-to-3D correspondence problem mentioned above.

6. Determining "true" shape of dents
Dents are hard to measure. A sharp dent and a rounded dent of the same surface area and depth will appear different. Measuring dents by mechanical means leads to subjective judgments, and it's also quite tedious to use depth gauges, rulers, etc., when you'll likely have to measures dozens if not more.

This is one of the harder engineering problems to solve for any defect detection project for manufacturing: how does one measure a defect and characterize it? if there is a standard for doing so, does the standard correlate well to something the inspection system measures? if the inspection system doesn't find a ding it "should have" found, who's to blame?

That said, if an inspection system works well enough for a sample of known defects, then users may eventually come to trust it, and the system itself becomes the standard for defining defect severity.

7. Extensive database of dents
Ideally you would have hundreds if not thousands of sample images of dents of different severities at different locations on vehicles from difficult manufacturers. If you're interested in finding dents caused by accidents during the assembly process, then it could take a long time to collect that kind of data. Dents caused during the assembly process are not common.

If you're only interested in finding dents caused by accidents or environmental damage, then that's a different matter. The types of dents will be different from those caused by accidental bumps inside an auto assembly plant.

8. Dealing with different paint colors
It's true that edge detectors can be reasonably robust at detecting edges in images of varying levels of contrast, but it can be quite disheartening to see what "varying levels of contrast" really means for different automotive paints and finishes. A light stripe that looks great on a shiny black car could be hardly detectable on a white car with old paint.

Most cameras have relatively limited dynamic range, so achieving good contrast for both black shiny surfaces and white dull surfaces is tricky. It's quite likely you'll have to automatically control lighting intensity. That's hard, too.


compression algorithm Steganography can hide one image inside another.

If we have a compression algorithm that allows us to fit the bits for two images into the size of one, can we think of a way to modify a run-length encoding algorithm to take two images and produce one, in such a way that the first is retrieved in a lossy way, and the second is retrieved?

Please explain this with full description of your thinking.

If you know a good reference, please let me know.

I wrote this program long ago and published, but removed it later, because could be an abuse. It is possible to write a plugin for the browser, that extracts hidden images and show them in the browser in place of originals, so some visitors will be seeing art gallery and others, who downloaded plug in, something else, may be illegal. I don't mean regular pornography, it is legal. It may be used by terrorists.
The way I did it was follows: I compressed one image and redistributed it in two least significant bits of each byte of pixel data. I used my own compression algorithm. In steganography the compressed image is called payload and where it is embedded is called carrier. My payload was half of the area of carrier. And I hid even multiple images within each other like Russian dolls. The carrier can be compressed and hidden within new carrier. However each payload was half of carrier area. I tested those program that claimed to determine steganographic content. Neither of them detected payload. However, there is a way to find that out. The carrier must be lossless, I used PNG. And that raised suspicion. Why image is posted in PNG format. Steganography is not illegal, I decided to remove my program from internet myself, there was no pressure.


Estimating the Impulse Response of the Room Using Sweep Signal Microphone Recorded Signal (Input & Output of a Convolution) I played this signal A (a 20Hz to 20000Hz sinusoidal sweep in 10 seconds) with a studio monitor speaker in a big church, and I recorded the result B with good microphones.

The result is very reverb-ish, that's exactly what I wanted to catch.

Now a software (such as Deconvolver but non open-source) can build an impulse response from A + B, that can be later used in a convolution reverb.

It works well. But I would like to learn how to do this myself via DSP / programming.

How can I use a sweep (signal A) + recorded output (signal B) to get an impulse reponse?

Edit: In other words, if a is the original sweep, b the recorded output, and h the impulse response, how to get h from

$$a * h = b$$

Is this formulation correct? is the solution $h=a^{-1} * b$, where $a^{-1}$ is the inverse of $a$ for the convolution? How to compute a convolution-inverse of a discrete signal?

this is the two-channel FFT method of spectrum analyzer:

$$ y[n] = h[n] \ \circledast \ x[n] $$

just make sure that the length of the FFT $N$ is at least as large as the length of sound $x[n]$ plus the expected length of the impulse response $h[n]$.  the length of sound $y[n]$ is also as long as the FFT.  you can round $N$ up to the nearest power of two.  just zero-pad everything to that length and then

$$H[k] = \frac{Y[k]}{X[k]} $$

is functionally true.

you might sometimes have to worry about division by zero, but if your driving signal $x[n]$ is sufficiently broad-banded (which a linear sweep or a maximum-length sequence is), then you don't have to worry too much about division by zero.

if you know $H[k]$, then you know $h[n]$.


2D to 3D view conversion I came across this demo by Fox sports:
https://youtu.be/N9SxFf5WCb0
Where they show a player's view in a 360 space. Very smooth and nice.
I was wondering how do they achieve it without any cameras attached to the player?
One approach that came into my mind was using a multi-camera setup and convert 2D images to 3D using a stitching mechanism as a post processin. Is it what they do?! What are other approaches that we can do it?

It seems to be at a stadium which uses 12 cameras and Replay-Technologies special signal processing called FreeD: https://www.youtube.com/watch?v=N1kgt2VDjdM
The image you show seems to use 2-3 cameras, but the technology can work with multiple camera combinations from the 12 in 360'.
It's a particularly stunning real-time effect, and there are also very many 2D to 3D apps available for example in after effects and on google scholar.
Here is a list of replay-technologies videos:
https://www.youtube.com/results?search_query=replay+technologies+freed
https://www.tvtechnology.com/news/bringing-the-matrix-to-the-nfl


What's the relationship between the parameters of this noise generator and the gaussian curve it produces in the frequency domain? This algorithm produces a sequence $y$ having a magnitude spectrum that resembles a normal distribution, peaking at frequency $\omega$:

$$b[n]=(1-\sigma)b[n-1]+\frac{\sigma}{\omega} C$$
$$p[n]=p[n-1]+\omega(1+b[n])$$
$$y[n]=cos(p)$$

where $b[n]$ is a leaky integrator having feedback $1-\sigma$ and taking white noise $C$ as input, $p[n]$ is a phase accumulator having frequency $\omega$ modulated by brownian noise $b[n]$, and y[n] is the final signal, cosine at phase $p[n]$.

Python function:

def signal(f, w, n):
    """
    Arguments
    f - peak frequency in radians/sample
    w - bandwidth of frequency content
    n - length
    Returns
    sig - list of samples

    """
    br = 0
    p = 0
    sig = [0] * n
    for i in range(n):
        # br - leaky integrator having rate proportional to bandwidth,
        # integrating white noise at magnitude inversely proportional
        # to frequency.

        br = (1 - w ) * br + w * random.uniform(-1, 1) / max([0.001, f])

        # p - phase increments at frequency f modulated by brownian noise.
        p = p + (1 + br) * f

        sig[i] = math.cos(p)

    return sig

To illustrate:

averaged magnitude spectrum at swept frequency from [0,π)

averaged magnitude spectrum with swept bandwidth from [0,0.5)

fourier transform of a sample sequence, hue indicates phase

wav

audio sample at 1760hz with $\sigma=0.005$

wav

averaged eight times

There are undesirable characteristics:

The amplitude of the peak frequency has an unknown relationship to $\sigma$, the second animation shows the sharp curved response
$\omega$ is not a perfect predictor of peak frequency, $\sigma$
allows noise to perturb the phase towards a range of
$[0\omega,2\omega]$ from FM

Short of using optimisation and curve-fitting, how can the frequency content be predicted (e.g. by the gaussian function) and controlled?

Any additional insight appreciated

This seems to work fine if everything is in rad/s. I think that the "problem" is with what you have implemented.

Here is a back of the envelope implementation of what you present:

Fs = 1000; % Sampling frequency in Hz
T = 4; % Total duration of the output signal
f = 125; % Centre frequency in Hz
s = 10; % Bandwidth in Hz

% Pre-allocate the output
y = zeros(1,T*Fs);

% Simple "buffers" for the previous sample in each filter
pb=0;
pp=0;

% For the duration of the signal...
for n = 0:(T*Fs)
    % Integrate b
    b = (1-2.*pi.*s./Fs).*pb + ((2.*pi.*s./Fs)./(2.*pi.*f./Fs)).*unifrnd(-1,1);
    % Integrate p
    p = pp + 2.*pi.*f./Fs.*(1+b);
    % Produce the output
    y(n+1) = cos(p);
    % Store the current values to the "buffers" to produce the n-1 indexing
    pb=b;
    pp=p;
end;

% Frequency vector for the FFT assuming nbins=length(y)
f_v = Fs.*(0:(length(y)-1))./length(y);
% Spectrum of the produced signal
semilogy(f_v, abs(fft(y)));grid on; xlabel('Frequency (Hz)');ylabel('Amplitude');

For f=125, s=0:

Clean as a whistle, right on the f

For f=125, s=15.625

For f=125, s=62.5

Notice that I provide the bandwidth at fractions of f. The bandwidth cannot be wider than f.

I think that this is closer to the discrete Gaussian as it is more akin to a power law for very narrow bandwidths.

The deviation is at:

br = (1 - w ) * br + w * random.uniform(-1, 1) / max([0.001, f])

Why .../max([0.001,f])? The uniform() already returns numbers between $[-1..1]$ why do we need the additional "normalisation" there?

Overall, this is like an oscillator whose phase becomes a random variable. No matter how you are going to achieve it, this oscillator's phase must accumulate with an average rate of $\Delta \phi = \frac{2 \pi f}{Fs}$. The deviation from this would be proportional to the width around $\Delta \phi$. When setting a bandwidth, you are trying to limit the "phase growth" between two limits (but again, on average). So, for the Gaussian, the $\sigma$ (the standard deviation) does not correspond to the point in the curve that the output is at half the maximum. This occurs at $\approx \pm 1.1775 \sigma$, which you can apply as a correction factor to your s to modulate where you want your amplitude to be positioned approximately. But that would still depend on a few factors (for example the length of the signal and the whole sequence of random values that determines constructive and destructive interference throughout the signal).

Since you are in Python, you can construct such a "probabilistic" oscillator by drawing samples directly from a Gaussian $\mathcal{N}(\mu_f, \sigma_f)$ where $\mu_f, \sigma_f$ have frequency units as above) to obtain a series of $\Delta \phi$s and then simply add them up as you along (think of it as cumsum). Of course, this is not going to work if the hardware this oscillator is to run on can only generate uniform random numbers.

Hope this helps.


Find step response for LTI system represented by the following impulse response? (i) $ h[n]=\left(\dfrac{1}{3}\right)^nu[n]+\left(-\dfrac{1}{3}\right)^nu[-n-1]$

here's how i tried to solve :

since we know, step response = running sum of impulse response
i.e, $$s[n]=\sum_{k=-\infty}^{n} h[k] $$
we can write

  $$ s[n]=\sum_{k=-\infty}^{n}\left(\dfrac{1}{3}\right)^k\left[u[k]+(-1)^ku[-k-1]\right]$$
  $$s[n]=\sum_{k=-\infty}^{-1}\left(-\dfrac{1}{3}\right)^k +\sum_{k=0}^{n}\left(\dfrac{1}{3}\right)^k$$
  $$s[n]=(-\infty .......+3^4-3^3+3^2-3)+\dfrac{3}{2}.\left(1-\dfrac{1}{3}^{n+1}\right)$$
  $$s[n]=\left(-3+3^2-3^3+3^4......-\infty\right)+\left[\dfrac{3}{2}-\dfrac{1}{2}\left(\dfrac{1}{3}\right)^n\right]$$and since sum within parethesis is diverging therefore,
  $$s[n]=\infty$$
  but, answer in the book is given as $$s[n]=\left(2^{n+1}-1\right)u[n]$$

above s[n] denotes step response

please tell me where i'm doing wrong(or if i'm doing conceptual mistake)
,and provide solution if possible ....thank you.

The given impulse response is unbounded for $n\rightarrow -\infty$, and, consequently, the corresponding step response is also unbounded for all $n$, as you've shown. The step response given in the answer corresponds to a totally different impulse response:

$$h[n]=2^nu[n]\tag{1}$$

which is also unbounded for $n\rightarrow\infty$, but since the corresponding system is causal, the step response assumes finite values for finite $n$.


8*8 block matrix in JPEG image compression? In standard JPEG format of a image, Discrete Cosine transform is used. But instead of applying the transform on whole image, we first divide the image in 8x8 block and apply transform on each of them. Thus during quantization we remove small coefficient in higher frequency. These steps are explained in details here.

But as we know that in N-samples DCT we can have only N frequecy coefficient. If we apply the DCT on the whole image we will able to get much higher frequency coefficient.

Can someone here highlight any disadvantage of applying DCT over the whole image at a time?

The lossy JPEG compression does not merely remove small coefficients in higher frequencies. It encodes them with a precision relative to a (relatively crude) visual perception model; most notably, horizontal and vertical frequencies are not quantized with the same precision. And as in many compression formats, it essentially assumes that the data is locally stationary.
If you apply the  DCT over the whole image, and  you quantify DCT coefficients, this quantization will effect the whole image. Imagine an image with a background checkerboard pattern, and a small zebra in foreground. With a whole-DCT compression, the zebra is likely to loose its stripes, because their energy is negligible with respect to that of the checkerboard. Even more, as JPEG applies DCT on chrominances as well, with down-sampling, color coefficient quantization is likely to produce false colors at places where they do not belong. With   block size larger than  about $8 \times 8$, several meaningful image gradients (parts of object edges) are more likely to happen simultaneously in the same patch.
Meanwhile, one can get slightly better results with $16\times 16$ blocks. That is for the perceptual part.
There is a computational part too. Small non overlapping blocks are easier to process, and require less memory, a still expensive part of electronic devices. Early JPEG 2000, applying wavelets on the whole image, have failed adoption despite better results, partly because of their memory footprint. Now, with the  advent of GPUs, processing blocks is becoming quite attractive. Last, a power of two ($8=2^3$) and clever tricks make the block DCT very efficient. In modern coders, standard often use different $2\times 4$, $4\times 4$... block sizes.
Last, but not least, it is not modular, as you need DCT implemetations for each image size.  Although it is possible with software,  chip makers are unlikely to like the hardware part.
So, applying DCT over the whole image:

is detrimental to  non-stationary images, and waster some local orientation,
can be costly in terms of memory, and hardware,
is not quite modular nor computationally efficient.

Thus being said, very stationary images would be better compressed with whole-size  DCT.


Effects of sampling rate approaching the Nyquist frequency? (not undersampling) I understand that you need to sample at a frequency greater than or equal to the Nyquist frequency to reconstruct a signal/aliasing etc. I read though that a "more conservative" approach is to sample at 2x-5x the Nyquist frequency (Rampil '98, on EEG processing). Assuming signal bandwidth doesn't extend beyond $F_N/2$ (has been low-pass filtered), are there effects/advantages of sampling at higher frequencies than $F_N$? That is, if I have a signal with no components above 100Hz, is there an advantage to sample at 1kHz vs at 200Hz?

The context is that I'm resampling various signals (EEG, EMG, joint angles) to a common time basis and need to pick the target frequency. Processing time is fast enough to not be a huge concern, but I do want to decrease it somewhat.

Thanks!

Sampling at 2*Fn is suitable only for signals that are infinitely long and perfectly bandlimited to below Fn, which conditions don't exist in actual signals being sampled or resampled (given non-perfect low-pass filtering), or in any non-infinite length signals.

A finite length rectangular window will itself introduce higher frequency content to the signal being sampled or recorded.

Therefore, there are major advantages to a sample rate higher than 2*Fn, in terms of not corrupting the signal, especially if any higher frequency (close to Fn) spectral content near the beginning or end of your recorded signal is important.

Added:  Reconstructing a signal represented by data at a sample rate well above twice the highest frequency can greatly reduce the computational (or hardware) requirements of the reconstruction filter or interpolator.


Using sinc as a filter I'm trying to use a cardinal sine as a lowpass filter for a cosine signal with a fundamental of 1k.

In frequency domain, what really happens is that I'm multiplying two impulses (centered at -1k and 1k) with a rectangular pulse of width equals 2 (convolution is multiplication in frecuency domain). The result would be a constant zero function. However using the above code in Matlab I'm getting again the sinc function as output. Any help would be appreciated (attached the full plot).

D = 5;
f0 = 1000;
fm = 2*f0;    % nyquist reestriction

t = (-D:1/fm:D);
x = cos(1000*2*pi*t);
h = sinc(t);

x_p = filter(h,1,x);
plot(x_p);

You messed with f0/fm relationship. Try this code:

D = 5;
f0 = 1e3;
fm = 40e3;
t = 0:1:2*fm*D-1;
x = cos(2*pi*t*f0/fm);
h = sinc(-D*pi:1/10:D*pi);
x_p = filter(h,1,x);
plot(x_p); plot(x); plot(h);
figure; psd(x,2^11,fm);
figure; psd(h,2^11,fm);
figure; psd(x_p,2^11,fm);


Real-time delay between signals I currently have a real-time noise suppression algorithm based on two inputs from different microphones. It has passed most tests that I have thrown at it, however when the signal I am trying to retrieve has a slight delay one mic relative to the other it works poorly.

I have done preprocessing delay estimation in the input signal through cross-correlation just for testing and it works, however in the final algorithm I can't preprocess the signal as it is supposed to be real-time.

I am taking sub windows of signals and applying my algorithm, what I tried to do is apply the cross-correlation to each sub window, it takes the FFT of the sub windows does the calculation and then the IFFT just to get the delay each time. I later have to take the FFT of the sub window with a different windowing function for my filtering algorithm. This much extra processing will drastically reduce the sub window size my algorithm can keep up with in real-time.

The question is, is there a faster way to do this delay estimation? Maybe something that works with each input sample to get a better estimation each time? Something like a control loop?

Thanks in advance!

I can't comment because I'm a new user... have a look to these questions:

Estimating the time offset (delay) between two audio signals in real-time

Generalized Cross Correlation

Time delay estimation of oscilloscope signals using cross correlation


Is there any research work on the possible data structures that can be used to represent an architectural design of a building? I am looking for a data structure, that can be used for description of a floor plan of a building. Its not about the format, I am asking about any other way to describe the floor plan of building, rather than in traditional map-based structure, preferably in some sort of tree or graph based structure, which can be derived from a text based description of the house.

DXF in computer aided design should already allow this representation. I would suggest you to check the format and adapt your data structures accordingly.

Check the data structures in LibDWG and LibDXF.


How to Create Image Editing Filters I am trying to make filters for days but now I just want to know the experience of others rather than go on to wrong way and finally decided to put my query over here with full detail rather then short question of every day . What is the best way of making image filters like Instagram or aviary or many other photo editing application , Can we make it using OpenCV or we need to make them by using operation on bitmaps like simple color operations of giving weights to them.

How they made the filters?
What is the logic behind it?
How we can perform the same operation?
Whats the best way of making them?

I read Quora - How Does Instagram Develop Their Filters?, but it didn't give the best way to operate with the tool for android . I don't even know that this question is right for this forum or not , but I just want to consult the algorithms , the best way to do it for android and simple the way from which I need to start ? I saw simple OpenCV - How to Implement Sepia Filter operation but its not easy to find that type of values for other filters like which we want to make like other colored shades in which only little effect can change the whole look of the image and the last thing is it the helpful if we make filters in Photoshop and then follow the steps in our program.
Thanks

The word "Filter" in Image Processing word relates to "Neighborhood" based operation on pixels.

Filters in the context of Instagram and other Image Processing applications are combination of point operations (Works on single pixel, adjusting it values independent of its adjacent pixels) and filters which alters the pixel value according to its adjacent pixels.

If I were you I'd start with learning about Point Operations such as Brightness, Contrast, Hue, Saturation, etc...
Once you build some functions which alters those move on to filters.

Then you could combine them into "Effect".

For instance, Vignetting is easy to do:

Create a rounded mask. White out side a circle and black around it.
Smooth the mask using a Gaussian Blur.
Adjust the brightness of pixel according to the values of this mask.


Summed-area table vs Integral image? I have been going through few research papers around face recognition and I came across two texture extraction algorithms :

Summed-area table by Crow et al. ("Summed-area tables for texture mapping", 1984)
Integral image by Viola et al. (e.g. "Robust Real-time Object Detection", 2002)

Apologizes if I'm wrong, I found that generating texture map in both these algorithms are same.

  i.e., Texture map $I(x,y)$ = Sum of all gray level intensity values
  who lies left and above $i(x,y)$ of the original image.

Please clarify me, what is the main difference between these two algorithms?

The Sum-area table (SAT) was introduced in computer graphics, and as far as I know, was popularized in computer vision and image processing under the name Integral image.

Apart from a larger genericity of the SAT to objects of dimensions different to that of images, the concepts are the same to me.