baseband and passband modulation I have a simulink model over 2.4GHz.

In the transmitter part, the modulation  is O-QPSK modulated baseband instead of passband modulation.

In Matlab Help this is the reason of using baseband modulator instead of passband one .

But I still can't understand how can linear channel  be  represented as baseband system as it can be represented as RF passband system.

I confused also between baseband and complex envelope.,

At complex baseband (meaning the signal is represented as I+jQ), the signal will have the same magnitude and phase versus time as the passband signal at any carrier given the system is linear.  For example, the signal 1+j0 has magnitude 1 and angle 0; we can also have a carrier at 1 GHz (or any frequency) with the envelope of magnitude 1 and phase 0 (relative to a reference at the same frequency). Likewise shifting the baseband by 45 degrees: .707 + j.707 results in magnitude 1 and angle 45 degrees is the same as shifting the 1 GHz signal by 45 degrees. Any linear operations on the two will be the same and it is significantly easier to model the baseband signal of interest versus keeping track of all those carrier cycles!  This is not very different from observing a spinning bike wheel with a strobe light, if we set the strobe light to be the nominal rpm of the spinning wheel, we can readily observe small changes in phase from nominal.

If the idea of a "complex baseband" is confusing, please also review this post here:

Frequency shifting of a quadrature mixed signal

This is also demonstrated in the following figures:

The first picture shows a complex baseband signal both in the time domain as a phasor moving randomly in magnitude and phase, but bandlimited such that it cannot move quickly. In the frequency domain we see a spectrum with bandlimited white noise (the spectral density is constant up to the bandwidth of the system).

Next we see the concept of frequency translation, where shown that multiplying by $e^{j\omega_o t}$ simply translates our spectrum to a new carrier $\omega_o$ but is otherwise unchanged.  The first plot below shows the exponential and it's Fourier transform, illuminating the frequency translation property since the transform is simply a frequency translated impulse: multiplying by the exponential (spinning phasor) in the time domain is convolution in the frequency domain, convolving with an impulse in frequency does not change any of the signal characteristics except for changing it's carrier frequency.

The next plot below shows the original signal we introduced, both at a carrier frequency and then "down-converted" to baseband by multiplying with $e^{-j\omega_o t}$, a minus sign in this case to move the signal from $\omega_o$ to baseband.  Notice specifically in this picture that the magnitude and phase shown versus time, on the IQ diagram, is identical for both: Both show the magnitude and phase relative to their respective carriers- in the plot on the left the carrier frequency is at $\omega_o$ and the plot on the right the "carrier" is at 0.  The strobe light analogy I gave earlier can also help with seeing this more intuitively.

Since these waveforms are identical under linear operations, I would usually choose to model and describe the waveform at baseband rather than at any actual carrier frequency, meaning in a simulation for example, I am not simulating an actual carrier waveform of which the signal of interest is modulating its amplitude and phase, but would do that directly at complex baseband with the same result.


How to Formulate a Constraint Which Ensures All Variables Have the Same Sign I'm trying to include a constraint in my problem (to be solved by any convex optimization solver). Let {a,b,c,d ...} be a finite set of continuous variables. How to formulate a constraint which ensure all of these variables, in the same time interval, take either "positive or zero" or "negative or zero" values.

sincerely,

Mohan

The solution from the blog you linked goes as following (Coordinating Variable Signs by Paul Rubin, Web Archive):

  Someone asked me today (or yesterday, depending on whose time zone you
  go by) how to force a group of variables in an optimization model to
  take the same sign (all nonpositive or all nonnegative). Assuming that
  all the variables are bounded, you just need one new binary variable
  and a few constraints.

  Assume that the variables in question are $ {x}_{1}, {x}_{2}, \ldots, {x}_{n} $ and that they are all bounded, say $ \forall \; 1 \leq i \leq n, \; {l}_{i} \leq {x}_{i} \leq {u}_{i} $.   If we are going to
  allow variables to be either positive or negative, then clearly we
  need $ {l}_{i} < 0 < {u}_{i} $.   We introduce a new binary variable $ y \in \left\{ 0, 1 \right\} $ and for each $ i $ the constraints
  become:

  $$ \forall \; 1 \leq i \leq n, \; {l}_{i} \left( 1 - y \right) \leq {x}_{i} \leq {u}_{i} y $$

  If $ y = 0 $ every original variable must be between its lower bound
  and 0, nonpositive. If $ y = 1 $ every original variable must be
  between 0 and its upper bound, nonnegative.

  Note that trying to enforce strictly positive or strictly negative
  rather than nonnegative or nonpositive is problematic, since
  optimization models abhor strict inequalities. The only work around I
  know is to change "strictly positive" to "greater than or equal to $
> \epsilon $" for some strictly positive $ \epsilon $, which creates
  holes in the domains of the variables (making values between 0 and $ \epsilon $ infeasible).

The problem with above approach is that it will make a Convex Optimization Problem into Non Convex Optimization Problem because of the multiplication between variables.

For many solvers I guess it will make no difference. Yet if you use Convex Optimization solver or even Disciplined Convex Programming (DCP) based solver (Like CVX) it won't work.

What I suggest doing, in this case, is take advantage that you have only 2 cases. One where all variables are Non Negative and the other when all of them are Non Positive.

Then all you need is to solve the problem with 2 different sets of constraints (Forming 2 different problems). When the constraints are:

$ \forall \; 1 \leq i \leq n, \; {l}_{i} \leq {x}_{i} \leq 0 $
$ \forall \; 1 \leq i \leq n, \; 0 \leq {x}_{i} \leq {u}_{i} $

Then chose the solution with the lowest objective value.

The main advantage is while the solution using Binary Variable (Which probably will do the same) require special MILP / MIQP solver the approach I suggested, for Linear Least Squares problem only require Non Negative Least Squares solver which is easier to find (Usually faster).

Example with Code

Let's say we want to solve the following problem:

$$\begin{aligned}
\arg \min_{x} \quad & \frac{1}{2} {\left\| A x - b \right\|}_{2}^{2} \\
\text{subject to} \quad & \operatorname{sign} \left( {x}_{i} \right) = \operatorname{sign} \left( {x}_{j} \right) && \forall i, j \in \mathcal{S}
\end{aligned}$$

Where $ \mathcal{S} = \left\{ i \mid {x}_{i} \neq 0 \right\} $, namely the set of all indices for which the element of $ x $ isn't zero.
So basically we say all elements should have the same sign or be zero.

Let's say we have a building block for solving Non Negative Least Squares. For instance one could use MATLAB's lsqnonneg().
To solve the problem above one could solve it one time for $ A $ and the other for $ - A $. Then take the argument which minimize the objective (If we take the argument of $ - A $ we need to negate it).

So the MATLAB function would look like:

function [ vX ] = SolveLsSameSign( mA, vB )

[vX1, resNorm1] = lsqnonneg(mA, vB);
[vX2, resNorm2] = lsqnonneg(-mA, vB);

if(resNorm1 < resNorm2)
    vX = vX1;
else
    vX = -vX2;
end

end

I created a simple script to compare results to CVX using the Mixed Integer Quadratic Programming (MIQP) which is given by:

cvx_solver('Mosek');
% cvx_solver('Gurobi');

cvx_begin()
    % cvx_precision('best');
    variable vX(numCols, 1);
    variable varY(1) binary;
    minimize( norm(mA * vX - vB, 2) );
    subject to
        vX >= (1 - varY) * lowerBound;
        vX <= varY * upperBound;
cvx_end

Pay attention that it requires MOSEK or Gurobi to be installed in order to solve this MIQP problem.

The full MATLAB code is available on my StackExchange Signal Processing Q52099 GitHub Repository.


High energy coefficients in 1d DCT for audio I've been reading a little bit about how DCT works in image compression using the 2D transform, however I don't know how this applies to audio using a 1D transform.

Suppose I have an audio segment, I mean, some samples of a mono audio vector.

Ai = [-0.0288,-0.0304,-0.0354,-0.0406,-0.0469,-0.0518,-0.0647,-0.0831 ...]

And then I calculate 1D DCT on that vector.

R = dct(Ai)

The result would be something like this:

0.2143   -0.3158   -0.6354   -0.3877    0.5503    0.2330   -0.1690    0.0203

How should I interpret this result? Are the high energy coefficient the ones with a greater value?(I'm still a noob on this). Is the first value the DC coefficient? Is there a special way to read the vector (zig-zag or something)?

Any help will be appreciated, btw I'm working with matlab.
Thanks

When we calculate 1D DCT of a signal (vector) we will get the DCT coefficients of same length.

The first value in this is the DC coefficient and low frequency components are at the begnning, as the index increases the frequency also increases. And hence the last value in DCT coeffs is the highest frequency (ie. fs/2).

Zig zag scanning is done in 2D DCT to read the coefficients in the increasing order of frequency. But in 1D DCT there is no need of zig zag scan as it contains only one row (or 1 column), it is already in ascending order of frequency (starting from Dc value and ends at maximum frequency ie. fs/2).


How to initialize observation Matrix in Kalman Filter when there is no clear relationship between measurement and state? I am try to use Linear Kalman to do time series prediction. I understand that I have to define a model process matrix which indicate how system state evolve, and a measurement matrix H which convert state variable to measurement space. However in my problem, there is no clear relationship between my measurement and state, so there is no way to just give 1 to some element in H to pick up the same physical meaning in state variable X, and 0 to other irrelevant component.
Under this situation, how do I design H matrix?
I have read a paper, it propose method below which use some data to train all the matrix, including H. But I don't understand how to deduce? can some one shed light on this?

In order to get a better understanding I used the method proposed in the paper to build the measurement function $y=h(x)$ using a Matrix $H$ and the state transition function $x_{n+1}=f(x_n)$ using a Matrix $A$. I did this for a practical example to get a feeling how it works.
To be clear in the paper you cited they call the measurements $z$ instead of $y$. So their measurement function is:

Their state transfer function is:

If the physical relationship between the measurement $z$ and the state $x$ is unknown one can use the formulas from the paper to calculate $H$ and $A$ depending on past states and measrements.
The example I tried is from Plane with constant speed. A radar system measures the distance to a flying plane that travels with constant speed. Starting position is 30 km away. Starting speed is 40 m/s. The state Vector has two values: position and speed of the airplane. So the inital state would be: $x_1=[30000, 40]$. Because the radar only measures the distance, the measurement equation simply is: $z_n = \begin{pmatrix} x_{1,n} \\ 0 \end{pmatrix}$. Thats because only the distance ($x_{n,1}$) gets measured but not the velocity ($x_{n,2}$).
The state transition equation is $\begin{pmatrix} x_{1,n+1} \\ x_{2,n+1}\end{pmatrix} = \begin{pmatrix} x_{1,n} + \Delta t \cdot x_{2,n} \\x_{2,n}\end{pmatrix}$. The velocity $x_{n,2}$ stays the same between step $n$ and $n+1$ but the distance changes by the multiplication of the time step $\Delta t$ and the speed during that time step $x_{2,n}$.
Using those equations and the measurement values provided in the example I get the following:

Next we try to achieve the same without writing the functions $f$ and $h$ explicitly but we try to use the formulas from the paper you cited.
For that also the recorded states $x_{n}$ and $x_{n+1}$ are needed. Luckily they are also provided in the example:

Using the Formulas for $H$ and $A$ and the recorded traces it is possible to calculate them.

$H$ in the case of our example becomes a $1\times2$ Matrix. Thats because $z$ is a $1\times10$ vector and $x_n$ is a $2\times10$ vector. As the latter is transposed $H$ ends up being $1\times2$.
Quite the same way $A$ ends up being a $2\times2$ matrix. $x_{n+1}$ and $x_n$ are both a $2\times10$ matrix. As the latter is transposed by the multiplication one ends up with a $2\times 2$ matrix.
Now we can rewrite the equations $f$ and $h$.
$$x_{n+1}=f(x_n): x_{n+1}=A\cdot x_n$$
and
$$z_n=h(x_n): z_n=H\cdot x_n$$
If I put in the recorded values I get $H=\begin{pmatrix} 0.956 & 35.38 \end{pmatrix}$ and $A=\begin{pmatrix} 1.00 & 4.98 \\ 0.00 & 1.00 \end{pmatrix}$.
With that I rerun the Kalman filter and I get the following:

This looks very similar to the former result which was derived using a physical modell. So in this case the method from the paper seems to work quite fine.
I did this using Matlab. So in case you are interested in the code:
close all

initialStateGuess = [30000 40];
n=length(initialStateGuess);
dt=5;

x_state=[30182  30351.4 30573.3 30769.5 31001.5 31176.4 31333.2 31529.4 31764.3 31952.9;38.2    36  40.2    39.7    43.1    39  35.2    37.2    42.1    39.9];
x_state2=[30373 30531.6 30774.3 30968.1 31216.8 31371.5 31509.2 31715.4 31974.8 32152.4; 38.2   36.0    40.2    39.7    43.1    39  35.2    37.2    42.1    39.9];
yMeas = [30110 30265 30740 30750 31135 31015 31180 31610 31960 31865];% Measurement Values
H=(yMeas*x_state2')/(x_state2*x_state2');
A=(x_state2*x_state')/(x_state*x_state');

% h=@(x)[x(1)];% Measurement function
% f=@(x)[x(1)+dt*x(2);x(2)];% State transition function

h=@(x)[H(1)*x(1)+H(2)*x(2)];
f=@(x)[A(1,1)*x(1)+A(1,2)*x(2);A(2,1)*x(1)+A(2,2)*x(2)];

ukf = unscentedKalmanFilter(...
    f,... % State transition function
    h,... % Measurement function
    initialStateGuess,...
    'HasAdditiveMeasurementNoise',true,...
    'HasAdditiveProcessNoise',true);

yTrue=linspace(30000,32000,10);
R = var(yMeas-yTrue);%0.2; % Variance of the measurement noise v[k]

Q=var(x_state-yTrue);
ukf.MeasurementNoise = R;
ukf.ProcessNoise = 30;

Nsteps = length(yMeas); % Number of time steps
xCorrectedUKF = zeros(Nsteps,n); % Corrected state estimates
PCorrected = zeros(Nsteps,n,n); % Corrected state estimation error covariances
e = zeros(Nsteps,1); % Residuals (or innovations)

for k=1:Nsteps
    % Let k denote the current time.
    %
    % Residuals (or innovations): Measured output - Predicted output
    e(k) = yMeas(k) - h(ukf.State); % ukf.State is x[k|k-1] at this point
    [xCorrectedUKF(k,:), PCorrected(k,:,:)] = correct(ukf,yMeas(k));
    predict(ukf);
end

plot([5:5:50],xCorrectedUKF(:,1))
hold on
plot([5:5:50],yMeas)
plot([5:5:50],yTrue)
legend('Kalman FIlter state','Measurement Values','True Value')
grid on


Autocorrelation of power spectrum Anyone have an idea of how I can implement autocorrelation of power spectrum of one image? I tried using:

 autocorrel = ifft( | fft(power spectrum) | ^ 2 );

but it's not working. Do you have any idea why or some other way?

The autocorrelation function is the inverse Fourier Transform
of the power spectrum. I have no idea where you got your formula
from or what your calculation is
meant to be doing. Are you sure that you are reading the
book/journal article/application note/website/Wikipedia or whatever
correctly?

In any case, if you want to compute an autocorrelation
as an ifft, you should be careful about making sure that
the power spectrum vector is set up right. If the
autocorrelation function that you compute as the iffy
does not have a peak at the origin (or as the first
entry in the vector for MATLABi folks) and have the
symmetry property $R[i] = R[-i]$ in real life and
the appropriately translated version in MATLAB, something
is awry.


Ouput SNR of an integrator Say I have an integrator with transfer function $\dfrac{\mathrm{F_s}}{s}$.

This integrator has an initial condition zero and integrates from $t=0$ to $t=\dfrac{1}{Fs}$.

If I apply a constant $V_{\mathrm{in}}$ together with a white noise source $e(t)$ ~ $N(0,\sigma^2)$ then I would have at the input $$\mathrm{SNR}=\dfrac{{\mathrm{Vin}}^2 }{\sigma^2 }$$

At the end of integration at $t=\dfrac{1}{Fs}$, the output power of the signal is $$P_{\mathrm{signal}} ={V_{\mathrm{in}} }^2 {F_s^2 \left(t\right)}^2 ={V_{\mathrm{in}} }^2$$
However, the output variance of the integrated noise (Brownian motion) is
$$P_{\mathrm{noise}} =\sigma^2 F_s^2 t=\sigma^2 F_s$$ resulting in SNR at the output Fs times lower than that at the input.

This calculation does not seem consistent with simulation.
Is there something wrong with my derivation?

I hope you do not mind, but I am going to change terminology a little bit. Since you mentioned simulations, suppose you are sampling the Gaussian white noise, at the input to the integrator, at a constant rate of $f_s$ samples per second. The point spacing, $Δt$, is $1/f_s$ seconds. Suppose $N$ independent consecutive samples are collected. Then the $N$ samples span a period of $τ_a$ seconds and this equals $N Δt$. The sample mean of the $N$ noisy inputs to the integrator obviously provides an estimate of the known $V_{in}$ input voltage. The sample variance of the $N$ noisy inputs to the integrator provides an estimate of the known $σ^2$ variance of the Gaussian white noise. The power SNR at the input is $(V_{in}/σ)^2$. The empirical estimate would be the square of the sample mean divided by the sample variance and bigger $N$ values are almost surely better.

The Gaussian white noise has zero mean and $σ^2$ variance, so the bilateral noise power spectral density, $η$, is $σ^2 Δt$. The ideal integrator has $τ_i$ integration time constant. (For an ideal op amp integrator, with just an input resistor, $R$, and feedback capacitor, $C$, $τ_i = RC$. The sign inversion is trivial and can be ignored.) Now suppose that the integration period, i.e., integration aperture or gate duration, is $τ_a$ seconds, exactly as above. Assume the integrator has initial condition zero, integration starts at $t = 0$ and ends at $t = τ_a$, so $τ_a$ is simply the reciprocal of the OP’s $F_s$. N.B. My $f_s$ is not the OP’s $F_s$.

Then the DC power gain of the integrator is $(τ_a/τ_i)^2$ and the noise equivalent bandwidth is $1/2τ_a$. Now consider two cases.

Case 1: Let $τ_i = τ_a$.  Then the signal power at the output of the integrator, at $t = τ_a$, is simply $(V_{in})^2$, because the DC power gain is unity. The noise power is $2η$ times unity DC power gain times $1/2τ_a$, which equals $η/τ_a$. But $η = σ^2 Δt$ and $τ_a = N Δt$. So the output power SNR is $(V_{in})^2 N/σ^2$. This is $N$ times larger than the input power SNR.

Case 2: Let $τ_i ≠ τ_a$.  Then the signal power at the output of the integrator, at $t = τ_a$, is $(V_{in} τ_a/τ_i)^2$, because the DC power gain is $(τ_a/τ_i)^2$. The noise power is $2η$ times $(τ_a/τ_i)^2$ times $1/2τ_a$, which equals $η τ_a /(τ_i)^2$. But $η = σ^2 Δt$ and $τ_a = N Δt$. So the output power SNR is $(V_{in} τ_a/τ_i)^2/[σ^2 Δt τ_a /(τ_i)^2]$, which again equals $(V_{in})^2 N/σ^2$. This is again $N$ times larger than the input power SNR.

In Case 1, the signal power is independent of $N$, while the noise power is inversely proportional to $N$, via $τ_a = N Δt$. So the power SNR is proportional to $N$. This is simply averaging independently sampled white noise. In Case 2, signal power is proportional to $N^2$, while noise power is proportional to $N$. So the power SNR is again proportional to $N$.

Relationship to matched filtering

The figure below shows that, if nothing was known about the input waveform for $t > τ_a$, then the input could be treated as though it was a rectangular pulse starting at $t = 0$ and ending at $t = τ_a$. The integration would effectively be gated integration, which would be a simple case of matched filtering: the impulse response of the (aperture synchronized) gated integrator is a rectangular pulse of duration $τ_a$ and reversing it changes nothing because the gated integrator's impulse response is symmetric. The figure below shows the calculation of the matched filter power SNR. The integration time constant, $τ_i$, has been set equal to unity because it affects both the pulse energy, and the noise, the same way, so it divides out.

I hope this helps!


What is limit on performance of edge-preserving filters for low signal-to-noise regimes? Let's say I have a signal, $x(t)$, defined as such,

$
x(t) =
\begin{cases}
 0 &\text{if} \,\,\, t < -\alpha/2 \\
 \frac 1\alpha t+\frac12 &\text{if} \,\,\, -\alpha/2 \leq t \leq \alpha/2 \\
 1 &\text{if} \,\,\, t > \alpha/2
\end{cases}
$

Thus, $x(t)$ is some signal that exhibits a jump from $0$ to $1$ over some timescale $\alpha$ around $t = 0$. For $\alpha = 0$, $x(t)$ has a discontinuous jump or an edge.

Now, let's add Gaussian Noise with $\sigma = 1$ to $x(t)$. Thus, $x(t)$ has a signal-to-noise of 1.

Now, say, I want to use some sort of edge-preserving smoothing filter on the noisy $x(t)$ to determine $\alpha$. It is obvious that any smoothing filter will wash away this edge a little and make it difficult to distinguish between a signal with $\alpha = 0$ and a signal with $\alpha = \alpha_{\mathrm{min}}$ for some minimum alpha. But there should be some minimum $\alpha$ that is distinguishable from a sharp edge.

For example, a Gaussian smoothing filter of some scale will make it impossible to distinguish between $\alpha$ less than that scale, but should make it possible to detect $\alpha$ larger than the scale of the filter.

What is the best possible performance a smoothing filter can have with regards to edge-preservation in this case? That is, what is the best possible smoothing filter that I can use for this low signal-to-noise regime that will be able to distinguish between $\alpha = 0$ and $\alpha = \alpha_{\mathrm{min}}$ for the smallest possible $\alpha_{\mathrm{min}}$.

Now, obviously, in practice we have discrete signals and $x(t)$ will be sampled at some sampling rate. So, hopefully, the answer to this question will address this. There's no real use to answering this question as if it was referring to a continuous signal.

This is Gibb's phenomenon (just not for the common edge case).

Hence, with a linear filter, the better you remove noise (which is, low-pass filtering), the more you remove the edges (which inherently have energy everywhere in the spectrum).

There's no minimum case – the noise bandwidth you allow is always the same as the signal bandwidth you allow, and thus, complementary, the amount of signal energy (and hence, edge precision) you cut off.

Among the things you can do is apply a clever signal model. Does your signal always look like your $x(t)$? Great, so you only need to estimate one parameter ($\alpha$) to fully describe it!


$1/n$ octave complex smoothing An excellent answer to this post explains how to do $1/n$ octave energy smoothing and mentions complex smoothing can be done as well, but it's tricky business because of phase wrapping.

How is phase wrapping solved in $1/n$ complex smoothing?

In smoothing, each data point is replaced by some kind of local average of surrounding data points. When the size of smoothing window is defined in terms of octave, e.g., 1/3th Oct., 1/6th Oct etc., it is called fractional-octave smoothing.

Let $Z[i]$ be the complex valued frequency spectrum, where $i$ is the frequency index (FFT bin) and $0 \le i \le N - 1; N$ is the length of frequency bins. Then smoothed spectrum is given by
$$
Z_s[i] =\sum\limits_{j=i-\Delta/2}^{i+\Delta/2} W[j] * Z[j]
$$
where $\Delta$ is the smoothing window and $W[j]$ is the smoothing coefficients. Here, $Z[j]$ is in complex form ($a + ib$) and $W[j]$ is real. So smoothed spectrum, $Z[i]_s$, contains both magnitude and phase information.

In general smoothing coefficients are normalized, such that, $\sum\limits_{j=0}^{n} W[j] = 1;$ where $n$ number of smoothing coefficients defined by $\Delta$.

There are different methods for choosing these smoothing coefficients. For example, Rectangular, Savitzky-Golay, Gaussian, Binomial, etc. All the methods are some kinds of moving weighted average. The idea is to give the most weight to the data points nearest the point of estimation and the lest weight to the data points that are furthest away.

In signal processing, the aim of smoothing is to enhance the S/N by reducing the noise as much as possible, but distorting the original shape as small as possible. There is a trade-off between reducing the noise and preserve the original shape.

You also have to handle the edging data points. There are several possible ways to deal with the boundary values, e.g., repeating the last value, zero padding, shortening the window size.


How would I find the function given the magnitude plot and the phase response? I'm wondering how I'd find the Fourier Transform X(jw) given the following information:

My understanding is that the expression for the continuous time fourier transform (CTFT) is magnitude(CTFT)exp(jw*phase), but I'm unsure how to use this expression here, since the phase is different for different values of w. I tried to do something piecewise:

from (-pi to 0):
X(jw)=rect(w/(pi))*exp(pi/2*j) = j*rect(w/pi)

from (0 to pi):
X(jw)=rect(w/(pi))*exp(-pi/2*j) = -j*rect(w/pi)

This doesn't agree with the answer key:
X(jw)=jrect((w+pi/2)/pi) - jrect((w-pi/2)/pi).

Would appreciate any thoughts on how to better approach this problem, since I'm very confused. Thanks!

***please note, here, j is used instead of i, the imaginary number (sqrt(-1))

There is more than one way to write the same thing. You can use $rect()$ functions or define it piece-wise, but in the end it is the same thing, since the $rect()$ is just a shorthand for a piece-wise definition.

If you do it in sections, you don't need to use $rect()$ function at all, you could simply write is as

$$X(j \omega) = \begin{cases}
 &  j \text{ if } -\pi < \omega < 0\\
 & -j \text{ if } 0 < \omega < \pi \\
 & 0  \text{ otherwise}
\end{cases}
$$

The answer just puts it together as two rectangle functions: one for negative frequencies and one for the positive ones.

Either way would be correct and I don't think one is "better" than the other unless there are some clear criteria what "better" means.


what is the the resultant, or, "final" frequency of the FFT of a chirped signal Let's say we have a time-transient with a chirped signal,
$$V(t) = \cos\left( 2 \pi f(t) t \right)$$
where $f(t)$ is the frequency of the signal and is time dependant. My question is how would one describe the final frequency, as produced by an FFT spectrum? My intuition is that it would be the time average such that
$$f^{\rm{FFT}}_{\rm{final}} = \frac{1}{T}\int^{T}_{0} f(t) \ {\rm{d}}t$$
where $T$ is the length of the time-transient signal used in the FFT. Strictly speaking I think would use a Gabor transform for non-stationary signals, but seeing as the spectrum analyser I am working with uses an FFT algorithm.

Depending on the chirp, if you seek a closed form solution, there may not be one. For empirical estimation, one can use non-stationary methods.
Take exponential with frequency ranging from 0.005 to .125 of sampling frequency. Applying synchrosqueezed CWT and plotting energy norm per row vs cumulative sum of rows' associated frequency:

The red line marks where energy contributions grow insignificant (and thus max freq is likely reached), which can be found with a flatline detection algorithm. It's heuristics and engineering.
Code below, can be improved.
import numpy as np
from ssqueezepy import ssq_cwt, TestSignals
from ssqueezepy.visuals import plot, imshow

def flatline_after_max(x, count=10, rel_dist=.1, max_frac=.99):
    # pick slightly less than max
    max_idx = np.where(x > max_frac*x.max())[0][0]
    ctr = 0
    for i in range(max_idx, len(x)):
        if abs(x[i] - x[i - 1]) / x[i] < rel_dist:
            ctr += 1
        else:
            ctr = 0
        if ctr == count:
            return i - (count - 1)

# define test signal ########################################################
N = 4096  # sampling rate (or number of samples)
fmin = .005
fmax = .125
x = TestSignals().echirp(N=4096, fmin=fmin*N, fmax=fmax*N)[0]

# synchrosqueeze & visualize
Tx, _, ssq_freqs, *_ = ssq_cwt(x, ('gmw', {'beta': 4}))
imshow(Tx, abs=1, yticks=ssq_freqs[::-1], ylabel="Fraction of fs",
       title="abs(SSQ_CWT)")

# compute row energies and find 'stopping point' #############################
Tx_energies = np.sum(np.abs(Tx[::-1])**2, axis=1)
Tx_energies_cs = np.cumsum(Tx_energies)
max_idx = flatline_after_max(Tx_energies_cs)
est_freq = ssq_freqs[max_idx]

# plot #######################################################################
plot(ssq_freqs, Tx_energies_cs, xlabel="Fraction of fs",
     title="cumsum(energy(Tx, axis=1)) | est_freq=%.3f" % est_freq,
     vlines=(est_freq, {'color': 'tab:red'}))


Using goertzel confuses me I am having some troubles understanding the algorithm.

My algorithm looks like this.

how come am i able to compute the DFT coefficient using this algoritm.
As far as i know
$$
X[k] = \sum_{i=0}^{N-1}x(n)W_N^k
$$

which using Goertzel is the same as $$X(k) = y_k(N) $$ and since $$y_k(n) = v_k(n) -W_n^kv_k(n-1)$$ and $$v_k(n) =2cos(\frac{2 \pi k}{N}v_k(n-1)-v_k(n-2)+x(n)$$ which means
$$X(k) = y_k(N) = v_k(N) -W_n^kv_k(N-1)$$ which isn't the same as $$X[k] = \sum_{i=0}^{N-1}x(n)W_N^k$$ ??

how come is Goertzel using $$|X(k)^2|$$ to compute the amplitude spectrum, where DFT uses $$|X(k)|$$ to compute the amplitude spectrum. ??

It looks to me as if you are asking two completely separate questions here.

First question: Indeed the two formulas aren't identical; it isn't obvious that the Goertzel algorithm computes what it does. You can find a sketchy proof based on the $z$-transform on the Wikipedia page about the Goertzel algorithm, or you can do a more direct proof by mathematical induction. Or you can just take it on trust -- for most practical purposes all that matters is that it does what it does, rather than why.

Second question:

With both the Goertzel algorithm and the complete DFT, |X|^2 is the power at a particular frequency. This is a sensible thing to compute because, e.g., when you do a DFT the sum of all the powers (|Fourier coefficients|^2) equals the total power in the signal (= sum of |signal values|^2, up to a normalizing constant factor).

If the frequency you're using is one of the ones processed by a DFT, the output of the Goertzel algorithm is (aside from "numerical" considerations like roundoff error) exactly the same as one bin of the DFT.


efficient Noise cancellation & suppression algorith to use in VoIP for low latency I am working on Noise suppression and cancellation in VoIP media engine. I am looking for noise cancellation & suppression algorithms and its implementation available with real-time low latency algorithms.

I had worked with webRTC noise cancellation implementation and it works well, but I want to explore more on this. I googled this and am not getting any concrete output. I need help from the community and experts to get better an idea of available noise cancellation and suppression algorithms. Below are findings:

Noise cancellation available implementation:

WebRTC
Speex

Noise cancellation available Algorithm:

Least Mean Square Algorithm

finding more on Github is giving lots of results, but needs advise from experts to the way to proceed for choosing Noise cancellation & suppression algorithm for VoIP with low latency.

You could use DSP algorithms based on deep learning techniques like this scenario:

You have clean voice Sig.1 .

You have known different types of noise Sig.2 .

Merge Sig.1 and Sig.2 then Sig.3 generated (noisy speech).

Feed it Sig.1 to deep learning based software also Sig.3 then make training the system.

Then obtained coefficients should be close enough to Sig.2

In this way you will have a system that could predict the noise type in real time and generated the anti noise signal which it opposite in phase (180) and add the microphone signal before send .


Matlab signal processing removing white noise from signal What is the way to address the removal of white noise from a signal?

Extra Details
Hi I also agree the details are as follow:                                                                                The signal is a acoustic or a sound wave. The the power associated with the frequency of the noise is distributed evenly across the whole frequency range (white noise) The actual signal characteristic is not known as its a inverse problem. But the actual frequency of the signal is a part of the total frequency. Power associated with the signal is higher than the noise can be seen in a spectogram.
Actuallty I also cant figure out how to extract the coefficients of the power/frequency  (db/Hz) from the spectogram. Can you suggest a way?
 see the graphs in the link

http://www.mathworks.com/help/signal/ref/spectrogram.html?refresh=true

Actual problem is to find the time of arrival of the actual signal

Without further details on the signal and the type of processing you want to apply, the two most basic approaches are:

signal enhancement via linear filtering (filter or filtfilt),
Wiener filtering, assuming a known stationary signal and noise spectra in an additive noise (matlab code).

If your signal is non-stationary, a time-frequency (spectrogram) or time-scale (wavelet) decompositions might help. In their domain, signal and noise are often better separated, and you can there apply Wiener filering, or alternative forms of thresholding, or more involved source separation techniques.

Wiener filtering in the wavelet domain are generated a lot of works, and codes, that you can test.


How is it possible plotting frequency vs. time in audio? Something which makes me confused is how it is possible to create plotting of frequency vs. time domain.

Isn't time domain is opposite of frequency domain and the only way is to plot amplitude vs. frequency or amplitude vs. time?

What is the main target of this plotting time vs. frequency of audio signal?

For traditional signal processing, time and frequency are dual variables, linked through Fourier transformations. This strong linkage yields a balance, sometimes called Heisenberg–Pauli–Weyl uncertainty inequalities. In everyday words, one could not be arbitrarily precise in time location and frequency determination. One cannot determine the frequency of a single sample.

However, in audio, meaningful information is usually carried by chunks of time, or collections of samples. So very often, people perform Fourier-like transformations over those consecutive or overlapping time chunks. They obtain information on how frequency contents evolve over time. A lot of specific methods, termed time-frequency analysis, have thus been developed, to provide  us with useful algorithmic tools to perform such analyses. They are legion: oversampled filter bank, spectrogram, Wigner-Ville, Choi-Williams, reassigned versions, S-transform, wavelets, etc.

In the tutorial: Time-Frequency Toolbox For Use with MATLAB (150 pages, no less), there are many comparisons. On a piecewise regular signal (sine then tone and linear chirp), they provide the following comparative plots that you can reproduce for your purpose:

One uses them  in audio for several reasons: to  analyze stationary parts, detect the onset of certains frequencies or harmonics, detect hearable information. Such techniques are at the very heart of mp3 audio compression for instance.


$\sin (t \omega)$ is not an Energy Signal, then how come its Fourier transform do exist? The following integral (perhaps fourier tranform of $\sin (t \omega)$ ) is not convergent:

$\int_{-\infty }^{\infty } e^{-i t \omega } \sin (t \omega ) \, dt$

As, $\sin (t \omega)$ is NOT an Energy Signal (but a Power Signal), then how come we get successful in finding the fourier transform of $\sin (t \omega)$ ?

You are right that such integrals are meaningless unless they are interpreted as distributions. And this is what we need to do, because - as you know - the Fourier transform of a sine function involves delta impulses. Let me try to make this a bit more intuitive:

The inverse Fourier transform of the delta function $\delta(\omega)$ (in the frequency domain) is given by

$$\frac{1}{2\pi}\int_{-\infty}^{\infty}\delta(\omega)e^{i\omega t}d\omega= \frac{1}{2\pi}e^{i0\cdot t}=\frac{1}{2\pi}$$

So we have the Fourier transform relation (time domain $\Longleftrightarrow$ frequency domain)

$$1\Longleftrightarrow 2\pi\delta(\omega)$$

Using the shifting property we obtain

$$e^{i\omega_0 t}\Longleftrightarrow 2\pi\delta(\omega-\omega_0)$$

And since

$$\sin(\omega_0t)=\frac{1}{2i}[e^{i\omega_0 t}-e^{-i\omega_0 t}]$$

we get for its Fourier transform

$$\sin(\omega_0t)\Longleftrightarrow \frac{\pi}{i}[\delta(\omega-\omega_0)-\delta(\omega+\omega_0)]$$


FFT for distance-limited sum of pairwise products I have a sequence of fixed length $N$ input arrays ($x_0$, $x_1$, ..., $x_{n-1}$), and my goal is to compute output arrays as such:

$$d_{00}[l] = \sum_{i = 0}^{N - l} x_0[i]x_0[i + l]$$
$$d_{01}[l] = \sum_{i = 0}^{N - l} x_0[i]x_1[i + l]$$
$$d_{02}[l] = \sum_{i = 0}^{N - l} x_0[i]x_2[i + l]$$
$$...$$
$$d_{10}[l] = \sum_{i = 0}^{N - 1} x_1[i]x_0[i + l]$$
$$...$$
With $l$ having a maximum value that scales linearly in $N$. That is to say, $d_{ij}$ will contain, in it's $k^{th}$ entry, the sum of all pairs in $x_i$ and $x_j$ that are distance $k$ away.

My understanding - Clearly, since each $d_{ij}$ output requires $l \propto N$ entries, each containing $N$ multiplies, the time complexity of producing each $d_{ij}$ output array is $O(N^2)$. These pairwise multiplies hidden behind a summation seem optimal for a Fourier Transform to improve the time complexity to $O(N \log N)$. However, I understand the standard Fourier Transform convolution does not allow for limiting to constant distance, instead allowing easy sum of pairwise indices in the form $\sum_{i=0}^{c} x[i] y[c - i]$.

How can I use this standard convolution, or any variation of the Fourier transform algorithm to make this problem computationally tractable?

Your sequences are all correlations between $x_i$ and $x_j$, where you assume that $x_i[n]=0$ for $n<0$ or $n\geq N$. Correlation is very related to convolution and it has a similar expression in frequency domain:

$$ c[n]=\sum_l x_1[l]x_2[l+n] \Leftrightarrow C[k]=X_1^*[k]X_2[k]$$

(capital letters are Fourier transforms of the small ones).

Hence, you can efficiently calculate the correlation using this theorem.

One caveat is that you are actually only looking for linear correlation, but using DFT you calculate circular correlation. So, the solution is to zero-pad the sequences. Have a look at the following Python code:

N = 5
x1 = np.random.randn(N)
x2 = np.random.randn(N)

x1z = np.hstack([x1, np.zeros(N)])
x2z = np.hstack([x2, np.zeros(N)])

# calculate according to your expression
d12 = []
for l in range(N):
    d12.append(sum(x1[i]*x2[i+l] for i in range(N-l)))

# calculate using zero-padding of the sequences with no
# dependence of the range of i
d12z = []
for l in range(N):
    d12z.append(sum(x1z[i]*x2z[i+l] for i in range(N)))

# calculate using circular correlation, using the zero-padded versions
d12c = []
for l in range(N):
    d12c.append(sum(x1z[i]*x2z[(i+l)%len(x2z)] for i in range(len(x1z))))
# check they all match
print ("x1:", x1)
print ("x2:", x2)
print ("d12_1:", d12)
print ("d12_2:", d12z)
print ("d12_3:", d12c)

# now, exploit the correlation theorem to calculate the same stuff:
C = (np.fft.ifft(np.fft.fft(x1z)*np.fft.fft(x2z).conj()).real)
# take only the first N values, as the others correspond to parts
# where the sequences do not fully overlap.
print ("FFT:", C[-np.arange(N)])

A sample output could then be:

x1: [ 1.11998979 -0.81341718  1.47117346  0.59699824 -0.9233909 ]
x2: [-0.31884542 -0.2214878  -0.02678967 -0.25943789  1.36355875]
d12_1: [-1.6303355800547898, 0.20609114482882795, 2.1870585244152387, -1.3997099032157734, 1.5271718773447336]
d12_2: [-1.6303355800547898, 0.20609114482882795, 2.1870585244152387, -1.3997099032157734, 1.5271718773447336]
d12_3: [-1.6303355800547898, 0.20609114482882795, 2.1870585244152387, -1.3997099032157734, 1.5271718773447336]
FFT: [-1.63033558  0.20609114  2.18705852 -1.3997099   1.52717188]


How to interpret an Ambiguity function I need to evaluate several different waveforms regarding their doppler tolerance, for use in a radar application. So I came across something called the Ambiguity function.

The following picture shows the Ambiguity function contour-plot of a Barker-Code with 7 Chips.

My problem is that I don't fully understand how to read it.

If I take the 0Hz Doppler-cut, I get the expected, ideal barker autocorrelation. Now consider a target which moves with a fixed velocity. In order to calculate the doppler frequency $f_d$ I can apply the following formula:

$$ f_d = \frac{\Delta v}{c} f_0$$

with $\Delta v$ being the (relative) velocity of the target, $c$ being the velocity of the medium (or the waves in the medium) and $f_0$ being the emitted frequency.

Now my problem is, that a barker code is not just built up by a single frequency which I could use as value for $f_0$. The barker code has a whole spectrum of frequencies (fft). This also means that I have to look at a bunch of doppler frequencies istead of looking just at one, right?....

So I cannot just take, for instance the 20kHz doppler cut, to see what correlation result I can expect for a target with a certain velocity.

What frequency of my signal-spectrum do I need to use to calculate the doppler frequency?

In my particular case I do a BPSK-Modulation of a 13x13 (nested) barker code, before sending it out and calculate the matched-filter response on the received signal. In order to get an idea of how the matched-filter output looks, I would like to take a "doppler cut" at $f_d = ..$. But should I take the carrier-frequency as $f_0$ ? Or should I take the maximal frequency that I have in my spectrum as $f_0$? ($f_0$ is needed to calculate $f_d$)

The ambiguity function is used in radar systems to get the distance and the relative speed of a moving object with respect to the transmitter. Is called ambiguity because it also tells about the ability to distinguish objects that are close between them and with a similar vector velocity.

The ideal ambiguity function is a Delta function in both domains. However there is no waveform with this ambiguity function. With your Barker-Code, you could have troubles to identify different objects that are at a distance +/- <10us, and whose Doppler is +/- <10 kHz. But you could distinguish two objects that are at 5us, one with 5 KHz and the other with 27 kHz.

In order to estimate the Doppler and the delay, you first compute the ambiguity function between the Barker-Code and the reflected signal. You will get a noisy version of the above figure but shifted in both domains, depending on the delay and Doppler of the target. If there are more targets, you will obtain overlapped ambiguity functions, each one shifted in a different place on the correlation plane.

You then just find the position of the maximum. From the obtained Doppler you can get the speed with the equation you wrote.


Why are some JPEG files of the same scene much larger than others? I'm using a Foscam FI8910W ip camera to view a static scene under constant lighting conditions.  When I pull back a frame grab, it's about 35 KB in size.  I can do this over and over and it's always around 35 KB, but fluctuates somewhat due to various noises inherent to electronic image capture.  This random fluctuation is only of the order of 1 KB at most.

About every 2500 frames, the image size of a frame is suddenly in the order of 70 KB.  There's no gradual creep upwards if you're thinking thermal noise as the camera warms up.  1 frame will be 70 KB(ish) and then they return to the 35 KB sized frames.

This has happened before with another run looking at a different scene.  The common file size was 39 KB then, and out of 10,000 frames, 4 were of the order of 77 KB.  The image size histogram looked like this:-

Before you ask, I've managed to save one of these frames and it looks exactly like all the others bar the expected noise fluctuation.  They have roughly the same number of unique colours at around 23,000.  So it's not a moth randomly landing on the lens for exactly 1 frame and then flying off.  For completeness, I've done another run of images and this is an example typical image (the reflection is the IR illuminator):-

This is the anomaly image:-

You can see there is no difference.  Excuse the hippo.  I'm fairly familiar with the JPEG algorithm, and I can't see how this might happen other than a coding error with Foscam's implementation.  But, could there be something inherently chaotic within some of the JPEG transformation functions (like the  discrete cosine transform or  quantization)?  Statistically, one would expect a normal distribution of file size and this is what I see around 39 KB.  Then there are several outliers at 77 KB.  So it doesn't appear stochastic.

The reason this is in CS and not in hardware is I'm asking could this be a programming code phenomenon relating to the JPEG encoding algorithm? Seems unlikely, but the anomalies are random and infrequent and there is no human interaction with the device.  Is JPEG encoding stable?

The reason you might be unfamiliar with this phenomenon is that as the images look the same, no one really looks at the file sizes.  The file size is of critical importance to me so I noticed.   How can this happen roughly every 2500 frames?

Supplemental:-

Posting these images is not going to work guys, as the imgur software re-samples the uploaded files.  So whilst I posted 37K and 73K files, imgur has re-sampled both to 35K.  This seems to be a Stack Exchange problem which seems ironic for a  site dealing with image processing, data compression and analysis.

This is my processing of the images.  It's the normalised difference between a normal image and the anomaly.  The image is as you would expect, with JPEG noise in the high frequency regions.    This is an RGB image even though it looks monochrome.  There are 8000 unique colours in the colour cube (representing the noise).

Supplemental 2:-

As requested, 4 normal frames and 2 abnormal frames can the downloaded from sample frames.  It is a different scene, but the anomalous behaviour still occurred, so this proves that it's consistent.

My guess is on autofocus or aperture briefly changing in a manner where the resulting image contains more high-frequency elements.

For example, if the focus moves from a smooth object to a textured object (like the smooth hippo to the cloth drape) in a manner where the latter's detail texture surfaces, JPEG will tend to take quite a larger size.

As someone else already stated: it seems like a good idea to check the image EXIF data for changes in core parameters like aperture and focus distance.  For such a striking difference in image size, it is quite likely that some fundamental parameter is different in the camera's opinion.


Frequency modulation implementation without trig functions https://dsp.stackexchange.com/a/1087 This answer suggests generating a sine wave using complex phasor multiplication. At each step the phasor is multiplied by another complex number (omega) to rotate it by the desired step. Then the imaginary part of the phasor is a sine wave. Can this be extended further to implement frequency modulation without having to use a trig function as part of updating the angular frequency? I've played with it a bit, but haven't figured out a clever way to remove the trig function when updating the frequency.

My initial idea: the phasor would be updated like so:

phasor = phasor * (omega * modulator)

The modulator would sweep back and forth between exp(j*freq_mod_max) and exp(j*freq_mod_min). I haven't figured out how the modulator would be implemented.

Side note: why would you want to do this? I did some quick and dirty performance tests and using this phasor multiplication method is quicker for generating a sine wave (about 3-4x quicker). Though calling the sine function directly is still of course quick enough so this is a somewhat academic curiosity at this point.

This is what an NCO (Numerically Controlled Oscillator) does. See this post for more details on the implementation and design considerations Numerically Controlled Oscillator (NCO) for phasor implementation?

By combining this with a D/A converter for an analog output we get what is commonly referred to as a Direct Digital Synthesizer (DDS).

With an NCO you can it you can easily do FM, PM and if you add an amplitude control at the output AM as well! To do FM specifically, you update the Frequency Control Word on each cycle according to the frequency vs time desired.

Since your question indicated a concern with efficiency; there has been a lot of work done in optimizing this approach to be highly efficient as well in terms of minimizing ROM storage. This includes more obvious approaches such as storing only 1/4 cycle of a sine wave, but also interpolation approaches are very common and other approaches include the Hutchinson Algorithm and Sunderland Algorithm.

The Hutchinson Algorithm uses a coarse ROM for a major phase step and fine ROM for a smaller increment and makes use of the small angle approximation:

$$\theta =\theta_F+\theta_C$$
$$sin(\theta) = sin(\theta_F)cos(\theta_C)+cos(\theta_F)sin(\theta_C)$$

for small anlge $\theta_F$:

$$sin(\theta) = \theta_Fcos(\theta_C)+ sin(\theta_C)$$

Thanks Hutchinson!

Also if you have time to spare in the implementation (meaning do not require a final update on each cycle), the CORDIC rotator is also quite interesting in terms of highly efficient approaches for implementing an NCO (or any other applications requiring phase rotation or trigonometric solutions).

References on the CORDIC:

Andraka, Ray, “A survey of CORDIC algorithms for FPGA based computers”, Proceedings of the 1998 ACM/SIGDA sixth international sympsosium on Field programmable gate arrays, pp191-200,  1998
Available at  http://www.adraka.com/files/crdcsrvy.pdf

The CORDIC FAQ at dspguru.com  http://dspguru.com/dsp/faqs/cordic

Maharatna, K, “Virtually scaling-free adaptive CORDIC rotator”  IEE Proceedings on Computers and Digital Techniques, 151, (6), 448-456.
Available at http://eprints.soton.ac.uk/263504/

Original foundation work on the CORDIC Algorithm:
Jack E. Volder, "The CORDIC Trigonometric Computing Technique," IRE Transactions on Electronic Computers, Vol. EC-8, pp.330-334, 1959.
John S. Walther, "A Unified Algorithm for Elementary Functions," Spring Joint Computer Conference Proceedings, Vol. 38, pp.379-385, 1971.
D. H. Daggett, "Decimal-Binary Conversion in CORDIC," IRE Transactions on electronic Computers, Vol. EC-8, No. 3, pp. 335-339, 1959.
M. A. Liccardo, "An Interconnect Processor with Emphasis on CORDIC Mode Operation," Masters Thesis EE Dept, University of California at Berkeley, September 1968.
J. E. Volder, "Binary Computation Algorithms for Coordinate Rotation and Function Generation," Convair Report, IAR-1 148, Aeroelectronics Group, June 1956.

en.wikipedia.org/wiki/CORDIC:  Matlab/Octave code for v=cordic(beta,n)


Band edge filter for raised cosine impulse I'm hoping to use an implementation of a frequency-locked loop for rough frequency synchronization in a PSK31 demodulator. The approach is to define a filter that is the derivative of the matched filter in the frequency domain. A couple slides from fred harris illustrate the idea:

I've found implementations for the typical case of a root-raised-cosine pulse shaping filter. However, PSK31 is a bit "special" in that it uses a raised cosine impulse. AKA, the Hann function:

$$ h(t) = {1 \over 2}\, (1+\cos(\pi t))\, \Pi(t/2) $$

Where $\Pi$ is the rectangle function.

How can I calculate the appropriate band-edge filter in this case?

It looks like the filter you want is indeed $-ith(t)$.  Here is some Octave code to get a visualization in the frequency domain:

t = [-1:0.01:1];
h = 0.5*(1+cos(pi*t));
hd = -i*t*0.5.*(1+cos(pi*t));
H=fftshift(fft(h,512));
HD=fftshift(fft(hd,512));
v = [-256:255];
plot(v, 20*log10(abs(H/512)), v, 20*log10(abs(HD/512)))

BTW, when plotting $h(t)$ in the time domain, your provided $h(t)$ appears to be normalized, spanning from $t = -1$ to $t = 1$.  So, I'm guessing your normalized symbol period is $T = 1$ with some ISI, or $T = 2$ with no ISI.


Counting the number of objects intersected by the image border in MATLAB I have an RGB image with various signs. My main goal is to count the signs that are in contact with the image borders.
Approach and problem
I started by loading the image [Fig. 1], then converted it to grayscale and applied a median filter to get rid of some noise [Fig. 2]. Then I binarized it with a threshold of 0.2, which resulted in Figure 3. At this time I got my binarized image, but the problem is that some parts that belong to the same sign are appearing in various regions, instead of only one. Now my goal is to merge the regions that belong to the same object , so then I could use bwlabel to count how many signs are in the image, and use imclearborder to get rid of the ones in the border, and use bwlabel again to get the difference between the two.
My approach was to use bwmorph, Dilate to dilate the objects and then try to fill them with imfill, holes. But the problem is that if I dilate them in a small amount[Fig. 4], the imfill doesn't seem to fill them, if I dilate them by a big amount[Fig 5] all the objects start to merge :(
Code
img=im2double(imread('image.png')); figure, imshow(img)
img_gray=rgb2gray(img); imshow(img_gray);
img_mediana=medfilt2(img_gray, [3 3]); figure, imshow(img_mediana);
img_bin=im2bw(img_mediana, 0.2); imshow(img_bin)
img_dilate=bwmorph(img_bin, 'Dilate', 10); imshow(img_dilate)
img_fill=imfill(img_dilate, 'Holes'); figure, imshow(img_fill)

Figures
Fig 1:
Fig 1 http://dl.dropbox.com/u/5272012/1.png
Fig 2:
fig 2 http://dl.dropbox.com/u/5272012/2.png
Fig 3:
fig 3 http://dl.dropbox.com/u/5272012/3.png
Fig 4:
fig 4 http://dl.dropbox.com/u/5272012/4.png
Fig 5:
fig 5 http://dl.dropbox.com/u/5272012/5.png

In my opinion, dilation erosion are very basic tools. You have a very strong information base and quite a decent input image to make such decisions.

Here is my view:

Given a reasonable success you have shown to go from Fig 1 to 3, you can identify and segment individual signs.
Assuming that you have known signs a prior, you can apply a decent fast algorithms to do pattern matching. In case if the exact patterns are not known- you can just identify the outer shape of the pattern.
Based on the classification, you can always define the centroid of each matched pattern and its respective width and height. If the centroid X,Y position is too close to the border - i.e. $centroid(x) < 0$ or $centroid(x) > imagewidth - shapewidth $ it is outside the edge, similarly you can apply for Y axis as well.
Given that you are only concerned about what falls on the edge - you should start with each edge only and start pattern matching there. Start matching the partial pattern/shape and if the partial pattern/shape does match that object IS being cut on the edge.

Here are some references that might help you formulate the problem well.

This paper is a very good to understand lot of basics about Signs/tokens you are dealing with.

Anil K. Jain and Aditya Vailaya Shape-Based Retrieval: A Case Study with Trademark Image Databases  Pattern recognition 1998, vol. 31, no9, pp. 1369-1390

There are many research elements which deals with partial or occluded shape/pattern matching.

Eli Saber,Yaowu Xu, A. Murat Tekalp Partial shape recognition by sub-matrix matching for partial matching guided image labeling Pattern Recognition 38 (2005) 1560 – 1573

Will expand this answer for more specific queries if you take this approach.


Poles and zeros of a transfer function What are the poles and zeros of this transfer function (in $z$):

$$H(z)=z+2+z^{-1}$$

and how would you approach the resolution of such problem?

Personally, I would write
$$H(z)=\displaystyle\frac{z^2+2z+1}{z}=\frac{(z+1)^2}{z}$$ then I would say that there is a zero (double) at $z=-1$ (numerator $= 0$) and a pole at $z=0$ (denominator $= 0$).

But wouldn't we also have a pole at $z=\infty$ since it makes $H(z)\to \infty$ which corresponds to the definition of a pole?
And in that case wouldn't that be contradictory since we know that $H(z)$ can be seen as a transfer function of a FIR filter which we know is always stable, so how can the filter be stable if it has a pole at $z=\infty$ (outside unit circle)?

I'm sure I have a basic misunderstanding about poles and zeros otherwise there shouldn't be any contradiction and hopefully someone can help me clarify this :-)

The "poles-inside-unit-circle" stability criterion only applies to causal systems. Your system is not causal because it uses one sample from the future owing to the $z$ term.

The general technique to check for stability involves looking at the regions of convergence (ROC) of $H(z)$. If the ROC includes the unit circle, then the system is stable.

See also this answer to the question: What is the relationship between poles and system stability?


Is Fourier series a sampled version of Fourier transform? I recently learned about dtft and how dft/dfs is the sampled version of dtft. I was wondering if Fourier series is also obtainable by sampling Fourier transform? I am a noob in the subject so sorry if this is stupid

There are 4 versions of Fourier transforms that are all close cousins.

It's all due to the basic property that "sampling in one domain corresponds to periodicity in the other domain". If a time signal is periodic, the frequency domain signal is discrete. If a frequency domain signal is periodic, the time domain signal is discrete.

So you have four variants

> `Time                       Frequency               Transform
Continuos & aperidoc          Continuos & aperidoc    Fourier Transform
Discrete  & aperidoc          Continuos &  peridoc    Discrete Time Fourier Transform
Continuos &  peridoc          Discrete  & aperidoc    Fourier Series
Discrete  &  peridoc          Discrete  &  peridoc    Discrete Fourier Transform

The flavors are different since a continuous signal requires integration and a discrete signal needs summation.

The naming conventions are non-intuitive: it would be easier if the DFT would have been called Discrete Fourier Series and the DTFT called DFT, but it is the way it is.

Which flavor you need to use depends on the nature of your signals and depending on what assumptions around the signals you can make, you can use one flavor to calculate a different one.


What are features in signal processing using a more tangible or less abstract description I'm learning about signal processing for music, and in my reading I keep finding the term "features". I find that they never really give a concrete example of what they mean by "features". They'll say stuff like features related to timbre, but they won't really go into further depth than that. I also took an image processing class once and found that I didn't really understand it there either. I understand what some literal features of a sound would be, like say distortion, but I have a feeling that this isn't really what they're talking about when they say features.

Example
I'm reading something on Dynamic Time Warping(an algorithm that seems similar to local/global alignment) and they say this

To compare two different features x, y ∈ F, one needs a local cost
measure, sometimes also referred to as a local distance measure, which
is defined to be a function c :F×F → R.

My understanding of Dynamic Time Warping is that two similar waves(∿), aligned to each other by measuring the distance between each subsection of one wave to each subsection of the other wave and then choosing the alignment that minimizes the total distance between all the wave sections. I'm doing this for soundwaves, so why are they describing it as comparing two different "features"?

Obviously it can't be described in a way that is 100% tangible, but I'm looking for a description of features that one would be able to draw out or make a picture of.

A feature is a number that describes one aspect of a signal.
Signals can be very complex, and the simplest analysis tools (like a time plot, a spectrum, or an energy measurement) don't tell you everything; in fact, for specific types of analyses, they almost don't tell you anything useful.
So, features are designed to describe very specific aspects of a signal.
To make an analogy, let's say you're an expert in ice cream and want to compare two brands. Just looking at a picture or the list of ingredients is not enough. You can taste them, but you can't accurately describe the taste to your colleagues (just like you can listen to sound but you can't describe exactly what it sounds like).
So, you come up with "features" about the ice cream: the density, the statistics of the nitrogen bubbles inside it, the sugar/fat ratio, the crystalinity, whatever. Some of the features will be useless, but some will turn out to be really good at objectively conveying something interesting about ice cream.


Order of operations when creating a spectrogram Signals noob here! I'm dipping my toes into your world by trying to recreate the Shazam algorithm, and I have a quick practical question:
When creating a spectrogram, is it better to take in a whole sound clip (say 5 seconds), and then chunk it into pieces for input to an FFT, or is it better to chunk it as the signal is being received and apply FFT as you go?
I'm asking from the perspective of ease of implementation, but if there are other considerations, I'd be happy to hear them.

To get a spectrogram using existing libraries, it is usually easier to do it for a whole auduio clip. For example in Python with librosa, it is a single call (assuming the audio is already loaded): librosa.stft.
A streaming implementation requires a bit more care. You will then need tohandling the window function, the FFT operation on each window, and the overlapped extraction of the windows. Not very complicated, but a bit harder.
The primary benefit would be lower-latency and lower memory requirements. If these are not key requirements, I would recommend the clip-wise solution.


Mahalanobis vs Normalization+Euclidean I'm using a set of features extracted from a signal for classifying the data window with KNN algorithm. Since the features have different value ranges, their influence on distance calculation is different when you use euclidean distance in KNN.

To equalize the influence of these features on classification:

I can normalize features with min-max normalization and then use euclidean distance.
I can use mahalanobis distance.

What are the advantages of these two approaches over eachother?

Both are reasonable approaches and it is foreseeable that either one could outperform the other empirically. The Euclidean distance assumes the data to be isotropically Gaussian, i.e. it will treat each feature equally. On the other hand, the Mahalanobis distance seeks to measure the correlation between variables and relaxes the assumption of the Euclidean distance, assuming instead an anisotropic Gaussian distribution.

If you know a priori that there is some kind of correlation between your features, then I would suggest using a Mahalanobis distance over Euclidean. Also, note that Z-score feature scaling can mitigate the usefulness of choosing a Mahalanobis distance over Euclidean (less true of min-max normalization though). The major drawback of the Mahalanobis distance is that it requires the inversion of the covariance matrix which can be computationally restrictive depending on the problem.


Time-Frequency Analysis of Big Data - Data Size Reduction: averaging the most appropriate method? Explanation:

I would like to analyse the data from an experiment, which investigates the performance of a mechanical component using sensors, that has generated 2000 CSV files. Each file contains 513 Rows x 1220411 Cols, and they are in spectrogram format (columns are time and rows are frequency):

| Time (s)| 0.0000 |0.000164|...
|:--------|--------|:-------|:-------
| 1.52kHz |  2747  |  350   |...
| 3.05kHz |  2996  |  420   |...
| 4.57kHz |  4078  |  300   |...
|   ...   |  ...   |  ...   |...

I have plotted my 3D chart of the first 100 rows using persp3D():

persp3D(x,y,z, theta=45, phi=5, xlab="Frequency (kHz)", ylab="Time (s)", axes=TRUE, expand=0.5, shade=0.2)

I would like to extract 1-4 column (1 sample of time for all frequencies) from each file to get a table of data with the total of 2000-8000 columns (for 2000 files with all the frequencies) and plot this to get the 3D plot of the experiment.

Question: best data reducing method

I would like to know what is the best method to make sure that those 1-4 columns represent the total data set in each file?
Is averaging simply a good method in this case? What are the alternatives?

So, first of all, CSV seems to me the least suitable format imaginable for this amount of data. It needs to be parsed, is memory hungry, and wastes precision, and isn't linearly addressable (ie. to get to the 99999. element, you need to parse the preceeding 99998 elements).

So I'd recommend keeping the structure from these files, but converting them to binary files, eg. numpy arrays, HDF5, whatever, but keep it binary, parse-free. Convert once and keep that data as what you load – it will reduce your file loading times immensely, and your software (whatever you're using) possibly doesn't even need to load all the file at once into memory, because when the data format has a fixed bit width per element (e.g. always 4-Byte integers and 8-Byte floating point numbers), it's trivial for the software to just load the part of the file into memory that you need. In fact, operating systems do that for you ("mmaped files").

Reducing your text file format to e.g. an array of single-precision floats or 32bit integers, your total data size would be a mere 8GB – and that fits multiple times in the RAM of every modern workstation, so I'd argue it's not really "big data" anymore! (Just to take the fear of handling so many numbers away)

I've ranted about this before, and people have answered that "CSV files are human-readable". I dare anyone who says so to read a 100k-column, 500 row CSV file and justify their impression of the data as being representable for the whole set.

  I would like to know what is the best method to make sure that those 1-4 columns represent the total data set in each file?

As always in engineering, there's no "best" without you defining a measure for good. If you want to highlight the differences among different frequencies, pick the column e.g. that has the highest sample variance. If you want to represent the whole measurement "fair", don't pick a single column, but average all columns. It's impossible to answer that question without being you – it's you that wants to show something. What you have is data – meaning is something that you humanly attribute (and, yes, you must be aware that the choice of data reduction method is your biased decision, and you should definitely communicate that).

  Is averaging simply a good method in this case?

It's the first thing I'd try to see if I can still see something interesting.
That's a pure coincidence – I might as well just pick a random one!

  What are the alternatives?

Impossible to tell, again – picking one column based on some metric, calculating variance, third statistical moment, postulate a certain functional shape would be handy and fit that shape to your data set … The sky is the limit, and what you want to find out is your only guideline.


Avoiding padding during FFT calculation I have a signal X(n) of length 24. But I want to find 1024 point FFT of this signal. I can do this by padding 1000 0s in x(n). But this requires huge memory which actually will be having only 24 valid data points.

Is there any way I can avoid this extra padding and yet be able to perform FFT. To be very specifica I am using CUFFT for finding FFT on the GPU. The padding 1000 0s means I have to allocate this extra memory on the GPU and also do copying operations from GPU to CPU and CPU to GPU etc.

Any hint?

From a programming point of view I don't see any possibility of avoiding the copy operation for 1024 values because the CUFFT library expects a pointer to an array containing NX (FFT size) values. If you don't initialize all of them, the result of the FFT operation will be random.

From a DSP point of view you could easily use a 32- or even 24-FFT without loosing information. Zeropadding in time domain corresponds to an interpolation in frequency domain, where the interpolated frequency bins can be different from the actual spectrum. You should check if you really need the frequency resolution of a 1024-FFT.


How to use IPP FIRGenBandStop method? Following this link:
http://software.intel.com/sites/products/documentation/doclib/ipp_sa/80/ipp_manual/GUID-5457048E-3C59-43E4-8A99-FB3873DCB769.htm
I have produced the following code:

        const int len = 512;
        Ipp64f pInputSignal[len];

        int tapslen = 27;
        int numIters = len;

        IppsFIRState_64f* pState;
        IppStatus status;

        Ipp64f* pOutputSignal = ippsMalloc_64f(len*sizeof(Ipp64f));
        Ipp64f* taps = ippsMalloc_64f(tapslen*sizeof(Ipp64f));
        Ipp64f* pDL = ippsMalloc_64f(tapslen*sizeof(Ipp64f));

        ippsZero_64f(pDL,tapslen);
        float sample_rate = 48000.0;//Hz //TODO DEFAULT_SAMPLE_RATE
        float frequency = 400.0;//Hz
        float amplitude = 2000.0;

        //Fill pInputSignal[t] = 2*pi*frequency*t, where t is time from t[i] = i/sample_rate
        fillBufferWithSinusoidData_64f(pInputSignal,len,sample_rate,frequency,amplitude);

        //computes tapsLen coefficients for bandstop  FIR filter
        float rFreqLow  = (frequency-52.0)/(sample_rate);
        float rFreqHigh = (frequency+51.0)/(sample_rate);

        status = ippsFIRGenBandstop_64f(
            rFreqLow,
            rFreqHigh,
            taps,
            tapslen,
            ippWinHamming,
            ippTrue);
        Assert::AreEqual((int)IppStatus::ippStsNoErr,(int)status);

        ////filter an input vector
        status = ippsFIRInitAlloc_64f(&pState, taps,tapslen, pDL);
        status = ippsFIR_64f(pInputSignal, pOutputSignal, numIters, pState);

In other words:
Input is a sinusoid at frequency 400Hz
Band stop between 348.0Hz to 451Hz

Expected output : something flat (as I would expect the sine signal would get filtered out).

Actual output : Sine wave moves 1/2 phase, magnitude is lowered a bit.

Any ideas about what went wrong here?

There are some issues with the code.
1. There seems to be a problem number of taps used for frequency. For filtering 400Hz sine wave at samplerate 48KHz I get better results if I use tapslen more than 559.
2. ippsFIRGenBandstop apparently is implemented as band pass + band stop + band pass. By doing this IPP tends to add some ripple at the higher end (around 44,000 and above). See here:

When using this filter it is recommended to use some kind of a FIR filter design tool like TFilter: http://t-filter.appspot.com/fir/index.html
Or matlab FIR filter design toolbox.


How are the scaling parameters included (extracted from) an essential matrix? An essential matrix relates corresponding points between two images assuming that a camera satisfies the pinhole camera model, and can be expressed as:
$$E=K'^TFK$$

where K is the calibration and F is the fundamental matrix.

There are many tools to compute the essential matrix, one of them using OpenCV's function findEssentialMat. You are typically returned a 3x3 matrix which contains the transformation (or the $$E=R[t]_x$$ camera pose version). Typically, the translation is directly extracted from the t part, and the rotation matrix has more than one solution.

However, how is scaling represented in the essential matrix?

In case it is, how can one extract it? In case it isn't, how come it isn't> And how should one include it to properly state the total motion between the images' points?

The Essential matrix is defined only up to scale, so you cannot extract scale from it. In other words, if you multiply $t$ and all the 3D world points in your scene by a constant factor, the essential matrix will be the same.

If you have to get the scale, then you need some additional information. Either you need to have an object of a known size in the scene, or you need some other sensor that tells you what the translation between the cameras is (e.g. GPS, odometer, or IMU). Here is an example where you get scale from a sphere of a known radius that happens to be in the scene.


Depth estimation from live video stream in matlab I want to estimate depth from live video stream with two calibrated webcams in matlab but has no results. In this post depth is estimated with two recorded video files. how can I estimate depth in live video stream?

Since you have two webcams you can reconstruct the video scene using the frame sequence. If your webcams happends to be stereo pair (with known extrinsic parameters aka rotation and translation) your task becomes pretty easy. On the other hand if webcam locations are unknown the problem gets a little trickier, but it can still be solved using matlab's computer vision toolbox.
In both of the above cases the scene reconstruction is accurate up to a scale factor. In order to calculate the depth in real units you need length knowlegde about an item in the scene. i.e the edge of a tracked square is 15cm.


Subpixel - what is it? I'm studying SuperResolution techniques, and I read often "subpixel"/"subpixel registration"/"subpixel accuracy". I don't find what exactly it means; Wikipedia says that is one of the three RGB channel in the pixel, but I think it isn't the right meaning in this context. Anyone can help explain what it is it or give me any tutorial?

Images are technically limited to the original pixel resolution. Objects in an image are spatially quantized at this resolution. However, an edge for instance may be localized "between" native pixels. Traditional tools like basic gradients often stick to that initial resolution. However, using additional information (multiple data, motion vectors, prior), one can infer a more accurate object location with a sub-pixel resolution, i.e. with an accuracy finer than that of the pixel width, at a fraction of the pixel size, with non-integer coordinates.

This can also apply to track/detect objects with a size smaller than that of a pixel. See for instance: What is a subpixel?

Using clever subpixel shifts, one can for instance reveal asymmetry not visible in the original image, as shown below:

Old and recent references:

Image Sequence Enhancement Using Sub-pixel Displacements, 1988
Multi-frame image super-resolution reconstruction using sparse co-occurrence prior and sub-pixel registration, 2013

Note: subpixels is also used in subpixel rendering on displays, and refers to colored pixel components (red, green, blue, plus  yellow, white  or cyan) that are displayed, or viewed, like a single pixel.


Library to differentiate people by their voice timbre Say there are 3 people in the room. Each of them have a unique voice timbre. I want to recognize people by their timbre and know who's speaking.

The task of taking a long contiguous audio recording and splitting it up in chunks in which only one speaker is speaking - without any prior knowledge about the voice characteristics of each speaker - is called "Speaker diarization". You can find links to research code on the wikipedia page.

If you have prior recordings of each voice, and would rather do classification, this is a slightly different problem (Speaker recognition or Speaker identification). Software tools for that are available here (note that general purposes speech recognition packages like Sphinx or HTK are flexible enough to be coaxed into doing that).


What is a 'notch' in image processing? I cannot find out what does "notch" means in image processing?
We have selective filters that one of the best of them is Notch filters. but what does the concept of "Notch"? why they are called "notch" filters?

In standard 1D signal processing, notch filters are filters that try to remove one specific frequency, and pass all others.

I say "try to" because this is not usually achievable, but it comes close.

Steve W. Smiths's DSP Guide has a nice plot (Figure 19.6b below) that shows the frequency response of such filters.

In image processing, it will be the same thing, but instead of a 1D filter it'll be a 2D filter probably comprised of two 1D notch filters --- one in each image direction.


How is a human whistle a pure sine wave? According to Wikipedia human whistling is a pure sine wave.

If that is true and if a human is whistling and produces a rhythm like any song, is he still producing a pure sine wave?

For example, if I am producing a tone of any music through my whistle would it still be a pure single harmonic sinusoidal wave?

This is only an approximation!

While the sound of whistling has a very prominent first harmonic, the other harmonics are still present too; there's a non-harmonic noise component; and the first harmonic is more a narrow bump than a sharp line. A better description of the sound of whistling would be white noise filtered through a band-pass filter with a very high Q, and sent through a slight waveshaper adding harmonics. With a slight modulation of amplitude and filter frequency.
A pure sine wave is a  stationary signal, it has constant frequency and constant amplitude. If you observe a whistling sound on a window larger than a few tenths of seconds, you will observe modulations in frequency (because of the melody and articulations between notes), and modulations in amplitude (as you breath and start/end notes).


Viterbi Algorithm - Selecting survivor paths which both contain the same Hamming distance To my (limited) knowledge, when determining which survivor path to choose to produce an error-corrected output of a binary Viterbi decoder, the path with the smallest Hamming distance is the most likely estimation of the transmitted signal.

I have come across an example where I am presented with two survivor paths, both with the same smallest Hamming distance. Is there a method of determining which path to then pick? Or is there no method of doing so?

Viterbi decoders are usually designed to take one of the two paths consistently versus the other in the case of a tie. This makes it a lot easier to do bit-exact non-regression testing of the implementation.

If the numeric precision of the input has enough resolution, the probability of a tie should be sufficiently small that it doesn't substantially impact performance to have this (presumably minuscule) bias. If there is not enough precision, then you will have non-trivial bias (and quantization noise), and such insufficient precision will  result in performance degradation, regardless of whether the bias is randomized or consistent.


Compensate for microphone? I have trained a machine learning model to recognize a specific audio event from mel spectrograms, and it works pretty well on a hold-out test dataset. I now want to run the model on a device and feed it with real audio but I have a problem that I don't now how to solve.

The problem is that when I'm testing my audio recognizer by playing audio clips it fails to recognize most of them, even if it can correctly recognize the audio in the original file. When I compare the spectrogram from the original audio file with the spectrogram computed by the recording of my device, I can see some differences:

Some of the differences has to do with the fact that the real time algorithm is doing the time-domain segmentation a little bit differently so that it "misses" the first few milliseconds of the event. But I have verified that this is not a big problem for the machine learning model. So the reason for why it fails seems to be because the loudspeaker and/or the microphone messes with the audio so that some frequencies are attenuated. The audio file has a 50kHz sample rate, which I think my high-end loudspeaker should be able to handle. My device samples the audio at 44.1kHz and the mel spectrogram uses frequencies up to 22.05kHz.

How can I make my feature extractor more robust? So that the spectrogram looks more similar irrespective of which microphone is used for recording or if I use a loudspeaker or if the audio event is actually happening in front of the microphone.

Some ideas: You could train the model using audio with randomly chosen speaker--mic system impulse responses applied. To get you started, there are microphone impulse response (IR) packs available, intended for music production purposes. Or you could perhaps generate synthetic impulse responses that cover the variability in real-world impulse responses. Or measure your own using various speaker--mic placements. The room will have an effect too; possibly you could use room IR packs, for reverbs.


monochrome (1-bit black-and-white) image conversion When converting scanned text to a 1-bit black-and-white image what are some filters that one can apply in the process to improve the result? Right now I am running into the problem where dithering errors make the image look horrible.

Update:
I think trying to undo dithering is a much harder problem. How can I convert the first image to a monochrome image? The default approach convert -monochrome img1 img2 is shown below. I've also tried a two-step approach: 1) decrease depth (color palette) and gamma 2) convert to a bilevel image (not shown). Other things I played with included imagemagick's ordered-dither (at various settings), but it wasn't as good as the two-step approach.

I would suggest for this case that you upsample and apply a slight blur and then a sharpen, then apply the threshold operation.  You won't get any more information from the pixel data, it's simply not there. But you'll get a smoother result out of the thresholding operation, and you won't need to dither. The end result is like a photocopier degradation.

Example:

Also, this is what it looks like when you just use a better diffusion dither algorithm ;)


what is the frequency response? I try to work out the simplest system, whose response is $H(z) = 1 + cz + c^2z^2 + c^3z^3 + \cdots = 1 / (1-cz)$

Now, Z-transform done quick (Fourier connection) and 4.5 Transfer Function, Poles and Zeros say that we just need to plug $z => e^{iw}$ to get the frequency response. I am getting single pole at $w = i \ln c$. It is complex. What is the frequency response? The second link says that real coefficients must generate complex poles in pairs. Where is the second pair?

Frequency response is $1/(1−ce^{iw})$, if you want to get amplitude information you just need to get the magnitude of this complex function and its phase will be the inverse of the phase of $1/(1−ce^{iw})$ or

$$-\arctan (-c\sin w/(1-c\cos w))$$

You can compute the magnitude of the transfer function by:

$$ 1 / sqrt( (-c\sin w)^2 + (1-c\cos w)^2 )$$


Why doesn't a sudden loud noise sound high pitched? The jumpiness or high change of a signal is due to the higher frequency components in the signal. So if I have a sound signal that increases suddenly in amplitude, why doesn't the signal sound very high pitched at that point?

For example, if I suddenly knock on a door or hit something, the pitch doesn't sound very high, but the magnitude of the sound increases from 0 to a number very quickly.

Pitch is a human psychoperceptual phenomena, which does not work the same way as a Fourier transform.

IIRC: Even if an FFT shows strong high frequency components in a short impulse or burst, the human ear-brain-system needs to detect roughly somewhere on the order of around 6 repetitions/periods/cycles of a periodic or pseudo-periodic acoustic waveform before communicating "I hear a pitch" to consciousness.  (depending somewhat on, but nowhere near linearly with, frequency). A single or duple period of even a pure sinewave won't do it.


Modelling room impulse by raytracing method I'm building a raytracing engine to get an approximation of the impulse response of a room, but the problem is that I don't really understand what to do after the ray tracing.

After I get all the reflections that collide with the receiver, I get all the rays with their acoustic energy and collisions time. But how do I turn this into an actual impulse response that i can hear using convolution. I was thinking in using a three band eq and turn the reflection energy into intensity, to set parameters for the eq, then send an impulse through the filter and do this for each ray; and finally sum all the impulse responses at their correct times. Please help, I'm not sure if the method I'm going to use is correct.

the impulse response is the sum of all of the different impulses received, beginning with the direct sound, then the primary reflections (only one reflection between source and receiver), then the secondary reflections (reflects offa two walls), etc.

the delay of the impulse is the delay of each path, which is the length of the path divided by the speed of sound.  each path has a path length and you get to apply the inverse-square law to that impulse delayed by that path length.  direct sound gets attenuated only by the inverse-square law.

all reflected sounds are attenuated by the inverse-square law and by each reflection.  you'll need to know something about the material of the wall at the reflection point to know how much attenuation.  a solid polished marble wall will have very little attenuation.


How to generate a noise like Pioneers white noise? I am making a music app and I want to put a noise generator in it. The problem is the usual white, pink noise is just too boring and not useful for the user. Pioneer mixers has an amazing noise generator. Check it here.
I want to generate something like this but obviously I have no idea how could such noise can be generated. I am new to DSP and all I need to know is what could be the magic behind this noise? Is it a High pass filtered white noise with some reverb? I need to get some guesses on the formula for this noise and I will do the research and implement it.

The "cool" stuff that you hear has almost nothing to do with the noise itself, but it's all in the post-processing: They use aggressive time variant low-pass and band-pass filtering, manual gain modulation and reverb/echo effects.
Any old pink noise will do just fine for this purpose.
"White noise" is an application specific misnomer. Almost no one uses actual white noise in audio, since it sounds really bad. Mostly it's pink noise or sometimes even brown noise.


How can I detect elements of GUI using opencv? Given a screenshot of an application, is there any way of finding GUI elements of that application using only opencv? (something like in the image)

I tried to play with blending the image and other effects from the coherence.py sample in the opencv installation folder. Did not succeed to get a viable result.

Has anyone did this? Do you have any idea what I might try?

Thanks

(I'm new to opencv and image processing so please don't avoid obvious explanations if possible) :)

OpenCV's MSER extractor (documented here) might be helpful — the bounding box of local MSER groups would pretty closely match the green rectangles in your mockup.


what does prototype mean in clustering? I am studying decision trees and machine learning. I have encountered with this word several time. Although, it seems that it not an important concept, I cannot understand this word in machine learning literature. For example:

  if the test does not significantly satisfy the criteria, then it creates a leaf and computes the prototype of the instances belonging to the leaf"
  .

what does prototype in this sentence mean?

According to Wikipedia:

  A prototype is an early sample, model, or release of a product built to test a concept or process or to act as a thing to be replicated or learned from.

But it does not make sense.

A prototype is an element of the data space that represents a group of elements. On the context of clustering (e.g. under a leaf), a cluster prototype serves to characterize the cluster, their elements. You can have a look at Cluster analysis: basic concepts and algorithms for instance, taken from Introduction to data mining.

In a more specific acceptation, it denotes an instance of clustering for combined numerical and categorical values, in the spirit of k-means, k-modes, k-medians, see for instance Into the world of clustering algorithms: k-means, k-modes and k-prototypes:

  For numerical and categorical data, another extension of these
  algorithms exists, basically combining k-means and k-modes. It is
  called k-prototypes.


How bandwidth is used in modulation In a typical Digital cable environment in the US, the spectrum is divided into 6MHz bands. When a tuner locks to a signal at say 800MHz, I would assume that the band for that signal is 797 - 803 MHz. The signal is QAM modulated. I would like to know why is that 6 MHz band required? If I can lock to a signal at 800MHz wont that be sufficient to reconstruct the original transmitted signal? From what I remember of my dsp course a few years back, I know that we need the band, but can't for the life of me remember why. Pointers to the theory will be appreciated

An ideal sinusoid does indeed have zero bandwidth, but it is not useful to transmit one because it also imparts zero information to the receiver.  Once the receiver knows what the sinusoid phase is it can perfectly predict the rest of the signal- thus, no more information is received.

For information to be transmitted there has to be an element of "not knowing".  For digital communications that generally means not knowing if the next bit will be a 1 or 0.  I will try to illustrate where the bandwidth comes from using BPSK, but the exact same principles apply to QAM.

The BPSK example above is ideal in the sense that it is very easy for the receiver to recover.  The receiver can sample the waveform anywhere within a symbol ("symbol" being the waveform representation of a digital bit) and get the full symbol power.  You can model the example BPSK waveform above as a bunch of individual square waves of varying widths and locations (in time) added together.  If you recall your Fourier theory, the Fourier transform of a square is the sinc function.  Thus, in the frequency domain, the total signal is a bunch of sinc functions of varying widths added together.  That makes the resulting waveform's frequency spectrum look a lot like a sinc function (note that this is not 100% true since there is some destructive interference that occurs, but for the sake of this discussion let's ignore that).

Sinc functions, unfortunately, aren't too great in terms of bandwidth.  The figure below shows the absolute value of a sinc function.

This doesn't look too bad until you look at the log scale of the energy, which is shown below.

As you can see from the plot above, at the location of the signal's potential nearest neighbor, the sinc function is only about 13 dB down.  The rolloff from there is quite slow, with it only being about 21 dB down at the fifth side lobe.  This introduces a lot of interference into neighboring signals, and when your receiver filters out the neighboring signals it is also filtering out a lot of your transmitted power, so that power is wasted.

Transmitters employ pulse shaping that rounds out the harsh transitions of the square waves and makes the BPSK waveform look more like a sinusoid.  The following is the same BPSK waveform but using raised-cosine pulse shaping.

And this is the energy vs. frequency view of the raised-cosine pulse.

The main lobe of the raised-cosine and the square wave are essentially the same.  The difference is in the rolloff of the side lobes.

So long story short (too late!) the bandwidth comes from the "random" transmission of ones and zeros.  The minimum bandwidth is determined by the symbol period (shorter period is faster, thus wider bandwidth), and the sidelobes outside of that main bandwidth are minimized via pulse shaping.


Anti-Aliasing and the Fourier Transform, Gonzalez Digital Image Processing In Gonzalez book Digital Image Processing, section 4.34 (third edition), he writes:

Unfortunately, except for some special cases mentioned blow, aliasing
is always present in sampled signals because, even if the original
sampled function is band-limited, infinite frequency components are
introduced the moment we limit the duration of the function, which we
always have to do in practice.
For example, suppose that we want to limit the duration of a
band-limited function $f(t)$ (i.e. a function whose Fourier transform
is non-zero only on a closed interval of range of frequencies), to an
interval say $[0, T]$. We can do this by multiplying $f(t)$ by the
function
$h(t)= 1 $ if $t \in [0,T]$, and is $0$ otherwise.
Then from the convolution theorem we know that the transform of this
product $h(t)f(t)$ is the convolution of the transforms of the
functions. Even if the transform of $f(t)$ is band-limited, convolving
it with $F(h(t))=H(\mu)$ will yield a result with frequency components
that are infinite.

This very last statement is what I am not sure about. If the Fourier transform of $f$ is band-limited, then outside of a closed interval, the transformed function will be $0$, and so I am not sure for which frequency components the convolution of the transforms will have infinite frequency. Any insights appreciated.

Note that in order to obtain the Fourier transform of the windowed time domain signal $f(t)h(t)$, you need to convolve the Fourier transforms of $f(t)$ and $h(t)$. We know that the Fourier transform of $f(t)$ is zero outside some interval (because $f(t)$ is band-limited). However, since $h(t)$ is time-limited, i.e., it is zero outside some time interval, we know that its Fourier transform cannot be band-limited. For this specific example of a rectangular window $h(t)$, we know that the corresponding Fourier transform is a sinc function. Convolving any function with a sinc function results in a function that extends from $-\infty$ to $+\infty$. Consequently, the Fourier transform of the windowed signal $f(t)h(t)$ also extends from $-\infty$ to $+\infty$. Consequently, the signal $f(t)h(t)$ cannot be band-limited.

In sum, a band-limited signal cannot be time-limited, and a time-limited signal cannot be band-limited. Note, however, that you cannot conclude that a signal is band-limited if it is not time-limited and vice versa. There are signal that are neither time-limited nor band-limited.