Implementing FIR/IIR "alike" filter in R I am new to signal processing. I am trying to simulate something similar to IIR/FIR filter with $k$ delays to imitate acoustic echo reflection. The difference equations for FIR and IIR respectively are as follows:

\begin{equation}
y(n) = x(n) + \sum_{D=1}^kA(n)x(n-D)+v(n)\;\;\;\;\; (1)
\end{equation}
\begin{equation}
y(n) = x(n) + Ay(n-D(n))+v(n)\;\;\;\;\;\;\;\;\;(2)
\end{equation}
where $D$ is a delay in samples, the coefficient $A(n)$ describes
the changing attenuation related from object reflection and
$v(n) ∼ N(0, 10^{−3})$ is the noise.

  Equation $(1)$ and $(2)$ can be found on section $VII$, second
  paragraph of this paper

How to can I implement this? I started by writing the following code in R.

Edit after the answer

install.packages("signal")
t <- seq(0, 1, len = 100)
x <- rnorm(100) + rnorm(length(t),0,0.001)
y <- filter(Arma(b = 0.1, a = 0.1), x)

Unfortunately, the above approach does not allow me to have $a=0$ or $b=0$.

Remark

Maybe following function is generating equation $(1)$:

fir1(39, 0.3)+rnorm(40,0,0.001)

The second equation perhaps is called flange IIR filter, where
the delay is not constant, but changing with time. This effect
imitates time stretching of the audio signal caused by moving
and changing objects in the room.

Response to the answer below

The ARMA equation

$$y(n)=\sum_{0}^Ma(m)x(n−m)+\sum_{k=1}^Kb(k)y(n−k)$$

i.e. AR part is as follows:

$$y(n)=\sum_{0}^Ma(m)x(n−m)$$

and MA part is as follows:

$$\sum_{k=1}^Kb(k)y(n−k)$$

which don't match with equations $(1)$ and $(2)$. In equation $(1)$ the coefficient $A$ depends on $n$. The coefficient update at each iteration.

Probably an answer

Parameters

t <- seq(1, 4000, by = 1)
x<- sin(2*pi*t*2.3)
A<- rnorm(4000)
v<- rnorm(4000,0,0.001)
k <- 20

Equation $(1)$

for(i in 1:4000){
  if (i>k){
    y[i]<- x[i]+ A[i]*sum(x[(i-(k-1)):i]) + v[i]
  }else{
    y[i]<- x[i]+ A[i]*sum(x[1:i])
  }
}

Equation $(2)$

for(i in 1:4000){
  if (i>k){
    y[i]<- x[i]+ A[i]*(y[i-i%%k]) + v[i]
  }else{
    y[i]<- x[i]+ A[i]*(y[1]) + v[i]
    }
}

Does the solution make sense? Or do I need to generate $x(n)$ using FIR and IIR filter?

The function filter() in the signal package does exactly what you are asking for. Note that, by default, it uses ARMA filters (a combination of AR and MA). Thanks to that, you can implement both FIR and IIR filters, because the difference equation for an ARMA model is:

$$y(n) = \sum_{m=0}^{M} a(m)x(n-m)+\sum_{k=1}^{K} b(k)y(n-k)$$

To achieve your FIR filter, note that it is enough to set all $b$ coefficients to $0$. For the IIR filter, you should set $a(0)=1$ and every other value of $a$ to $0$.

Another maybe obvious yet important observation is that you can add the noise $v(n)$ at the end of the script, because the noise is white (at least that's what I understand from your question). Thus, there is no correlation between past and present samples of the noise, letting you adding it after applying the filter.


Evaluate phase of $X(\omega)$ without computing $X(\omega)$ $$x(n) = \{ -1, 0, 1, 2, 1, 0, 1, 2, 1, 0, -1 \}$$
Let $X(\omega)$ be the DTFT of $x(n)$. I need to find the phase of $X (\omega)$ without computing $X(\omega)$. I notice that $x(n)$ can be a type I linear phase filter, but not sure if that is used in solving the problem.

hint: x(n) has symmetric coefficients which as a filter is a type 1 linear phase filter as the OP has stated. Such a filter will have a phase slope that is dependent on the number of taps only.


Would the locations of SIFT features generally agree with features detected by Shi-Tomasi method? I am having a comparison of different feature detection methods. So if we compare the SIFT (Scale-invariant feature transform) and Shi-Tomasi method, will their feature locations agree and why? Are they comparable? Thanks in advance.

Probably not. The SIFT detector finds centers of blob-like features. Shi-Tomasi detector finds corners. Furthermore, SIFT detector operates at multiple scales, while the classic Shi-Tomasi does not.


How can $a_n$ for the trigonometric Fourier series be $0$ for a non-odd square signal? Can $a_n$ for the trigonometric Fourier series be $0$ for a non-odd wave? Since at the moment the square wave is not even nor odd. Square wave is given by:
$$
s_2(t)=\begin{cases}
 -4, & 0\leq t\leq 2.5\\
6,& 2.5\leq t<5
\end{cases}
$$

The signal looks like this:

When I find $a_n$ it gives me $0$ but when you get a value of zero, doesn't this mean that the function is odd?
OR DOES THE ZERO ALSO IMPLY THAT IT IS NOT AN EVEN and clearly the wave isnt even?

EDIT:

If this wave is odd, how is it odd because I thought it had to be $s(t)=-s(-t)$ so for e.g. $6$ doesn't = $-4$?

Indeed, if the cosine coefficients ($a_n$) of the Fourier series of a function $s(t)$ are all zero, then the function is the sum of sine terms, which are all odd, implying that the function $s(t)$ is odd.

In your case, however, not all $a_n$ are zero; in particular, $a_0=2$. This means that the wave has a "DC" or mean value of $a_0/2=1$. Since $a_n=0$ for $n>0$, then the function $$u(t)=s(t)-a_0/2$$ is indeed odd, and now $a_n=0$ for all $n$.


How to modify an LTI differential equation to be acausal or anti-causal? I'm trying to wrap my head around causality in LTI-systems. Considering continuous time only, I'm happy with the fact that the system is causal iff the impulse response function $h(t)=0$ for $t<0$.

What I don't understand is how this is embodied in an LTI differential equation. Specifically, what makes the positive and negative time directions fundamentally different? (a point which is especially clear if considering e.g. length as the independent variable, rather than time)

I'm hoping to shed some light on this by considering the LTI equation

$ \frac{dy}{dt} = u(t)$

For this, or a more complex LTI system, I could set the RHS to $\delta(t)$ and then use Laplace transforms to solve for $y(t)=h(t)$. In this specific case $h(t)$ is the Heavyside step-function, i.e. $h(t)=0$ for $t<0$ and $h(t)=1$ for $t>0$. Thus this is a causal system.

If I wanted to modify the above differential equation to make the system either acausal, or anti-causal, what would I have to do?

It's important to realize that generally the differential equation (DE) alone doesn't tell us anything about causality. You claim that the system given in your question is causal. However, an anti-causal system with impulse response $h(t)=-u(-t)$ (where $u(t)$ is the step function) also satisfies the very same DE! So the given DE actually describes two systems, a causal one and an anti-causal one.

What you need to enforce causality (and linearity and time-invariance) is the auxiliary condition that the system is initially at rest. "Initially" means the moment in time when the input signal becomes non-zero. If we denote this moment by $t_0$, these conditions are

$$y(t_0)=y'(t_0)=y''(t_0)=\ldots =0\tag{1}$$

where $y'(t)$ denotes the first derivative of the output signal, etc. Solving the DE with the initial rest conditions $(1)$ for $t>t_0$ gives you the output signal of the causal system described by the DE.

A similar thing is true when using the Laplace transform to compute the transfer function of an LTI system. The expression for the Laplace transform of the transfer function does not uniquely identify the corresponding system. You also have to specify the region of convergence (ROC). In your example, the system's transfer function is

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

Depending on the ROC, this function corresponds to two different time domain functions:

$$h(t)=u(t),\quad \text{Re}(s)>0$$
and
$$h(t)=-u(-t),\quad \text{Re}(s)<0$$

You need to choose which solution makes sense in your case.


Increasing Image Resolution I know of some oscilloscopes (DSA8300) that repeatedly sample at a few hundred kS/s to reconstruct a few GHz signal. I was wondering if this could be extended to 2D signals (photographs). Can I take a series (say 4) of still pictures using a commercial 16MP camera to finally reconstruct a 32MP image? Will doing this remove aliasing I have from each image?

Should such a thing be attempted from a single image it would obviously not work as no new information is being introduced. If all the pictures taken are absolutely identical, will I still be at the same point as having one image? So are variations essential? Is CCD / CMOS noise enough variation for such a thing to work?

Is there a name for such a technique or algorithm? What should I look for?

One word for that technique is superresolution.

Robert Gawron has a blog post here and the Python implementation here.

Usually, this technique relies on each image being slightly offset from the others.  The only gain you'd get from not moving between shots would be to reduce the noise level.


FFT of square wave - what does output represent? I am really new to FFT and signal processing. I am doing an analysis of square waves with FFT and I am trying to understand why the FFT output on the frequency domain has a downward slope for square waves. Specifically:

Why in FFT output for a square wave, bin 0 has the highest amplitude and bin 1 has the second highest amplitudes no matter what the original frequency of the input signal in the time series domain representation is.
I read that bin 0 has to do with the ADC and that is the sum of the whole signal and that is why it is 0Hz. What does that mean?
What does amplitude represent in FFT output like below?

Simply put, a square wave can be decomposed into a weighted summation of series of sine waves, which is exactly what the Fourier transform does. Here are two animations that describe the synthesis of a square wave by accumulating sine waves.

The amplitude of FFT results are the weights of each sine waves. As seen in Fig.1, the first order harmonic has the largest weights, that's why the first bin has the highest amplitude, and the following bins have decreasing amplitudes.
An ideal square wave with an amplitude of 1 can be represented as Fourier expansion:
$$
x(t) = \mathrm{sgn}(\sin\omega t) = \frac{4}{\pi} \Big( \sin\omega t + \frac{1}{3} \sin 3\omega t + \frac{1}{5} \sin 5\omega t + \ldots \Big)
$$
The above equation shows that the ideal square wave contains only components of odd harmonics and their amplitudes are decreasing. However, in real world square waves contain all integer harmonics.
Bin 0 represents DC (direct current, i.e., 0 Hz). Since your square wave has a zero mean, the DC components is zero as well. Try to add a constant value to the original square wave, then take FFT and see what happens.


Detect signals degree of imbalance I have got this question: let's suppose that I have three signals, representative of the same variable, but taken on three different component of a same device; in my case, three load signals, representative of the out-of-plane bending moment Mx of each of the blade of a wind turbine.

Is it any possible to effectively and efficiently detect the degree of imbalance or asymmetry of such signals?

With asymmetry and imbalance I mean a single parameter/variable which gives information about how the three signals differ from each other.

I would then set-up a threshold value to this degree of difference, in order to further trigger a certain process when the threshold has been overcome.

The first idea I came up with is multiple correlation, but perhaps there exist tons of better methods.

EDIT

Each signal coming from each blade;
Blades placed with 120 deg phase shift;
Azimuth angle vs. Time available.
Regarding balance vs. unbalanced, I only know that unbalanced condition leads to higher loads on blades. In general, the three signals will always be unbalanced: what I would like to ascertain is the degree of imbalance; is it over a defined allowed limit or not?

One thing that can be useful in condition monitoring of rotating machines when rotation angle is important is to look at using circular or directional statistics.

Suppose your measurements are:
$$
m_1(n), m_2(n), m_3(n)
$$
for the three moments $m_i(n)$ for $i=1,2,3$.

Calculate the instantaneous directional mean of the moments:
$$
m_{\rm tot}(n) = \left ( m_1(n) + m_2(n) e^{j2\pi/3} + m_3(n) e^{j4\pi/3} \right)\\
\bar{\phi}(n) = \arg \left(  m_{\rm tot}(n) \right)
$$
assuming that the first blade is your reference.

This direction, $\bar{\phi}$, is the direction of the average moment.  If it's always aligned (or "anti-aligned") with one particular blade of your turbine, then that can indicate your "problem blade".

You can use the magnitude
$$
\left| m_{\rm tot} \right|
$$
as a measure of how "out of balance" the turbine is, as this should be (close to) zero most of the time.

If you have the azimuth angle, $\theta$, then you can get fancier. It may be that you need to angle-align (rather than time-align) your measurements to pick the faulty blade.


Discrete Fourier transform of even function I have a cosine function - Hence it is even. Considering only the real parts of the DFT, on performing DFT, I am getting something like this. Could anybody tell me where I am going wrong:

DFT:

Please let me where I am going wrong with the implementation

The real question is: what did you expect to see and why? I give you a few ideas and hints to help you understand the problem:

if you use the time domain data as shown in the top figure, then the first (left-most) point is interpreted as time index $n=0$ (there goes the cosine ...). Just look at the formula for the DFT and you'll see why.
Second, the implicit periodic continuation of the time domain data in the DFT means that you don't have a pure sinusoid (just append the first data point after the last data point and so on, and you'll see the problem). What you probably wanted is an integer number of periods in your time window (= DFT length), because then the periodic continuation will not change your nice sinusoid.


Histogram and binary image in opencv? I would like to compute how many zeros and ones (e.i. 255) in a binary image. The following code to generate a binary image by using Matlab.

I = imread('images.jpg');
level = graythresh(I);
BW = im2bw(I, level);
imshow(BW)
imwrite(BW, 'img.jpg');

This is my C++ code.

#include <opencv2\core\core.hpp>
#include <opencv2\highgui\highgui.hpp>
#include <opencv\cv.h>
#include <map>
#include <iomanip>

std::map<int, int> computeHistogram(const cv::Mat& image)
{
    std::map<int, int> histogram;

    for ( int row = 0; row < image.rows; ++row)
    {
        for ( int col = 0; col < image.cols; ++col)
        {
          ++histogram[(int)image.at<uchar>(row, col)];

        }
    }

    return histogram;
}

void printHistogram(const std::map<int, int>& histogram)
{
    std::map<int, int>::const_iterator histogram_iter;
    std::cout << "\n------------------\n";
    for( histogram_iter = histogram.begin(); histogram_iter != histogram.end(); ++histogram_iter)
    {
      std::cout << std::setw(5) << histogram_iter->first <<  " : " << histogram_iter->second << "\n";
    }
    std::cout << "------------------\n";
}

int main(int argc, char **argv)
{
    cv::Mat img = cv::imread("img.jpg", CV_BGR2GRAY);
    printHistogram(computeHistogram(img));
    return 0;
}

This is what I got

------------------
    0 : 11451
    1 : 590
    2 : 602
    3 : 498
    4 : 428
    5 : 376
    6 : 387
    7 : 349
    8 : 314
    9 : 278
   10 : 229
   11 : 206
   12 : 212
   13 : 159
   14 : 142
   15 : 112
   16 : 106
   17 : 96
   18 : 86
   19 : 54
   20 : 44
   21 : 37
   22 : 26
   23 : 26
   24 : 14
   25 : 14
   26 : 10
   27 : 4
   28 : 7
   29 : 8
   30 : 2
   32 : 7
   33 : 2
   34 : 3
   36 : 2
   37 : 4
   38 : 1
  211 : 1
  216 : 1
  217 : 1
  219 : 4
  220 : 2
  221 : 2
  223 : 3
  224 : 5
  225 : 2
  226 : 4
  227 : 4
  228 : 5
  229 : 8
  230 : 15
  231 : 19
  232 : 22
  233 : 25
  234 : 36
  235 : 21
  236 : 54
  237 : 54
  238 : 64
  239 : 93
  240 : 113
  241 : 111
  242 : 127
  243 : 163
  244 : 193
  245 : 182
  246 : 202
  247 : 218
  248 : 263
  249 : 289
  250 : 291
  251 : 301
  252 : 403
  253 : 458
  254 : 469
  255 : 29211
------------------

Is this right? What I'm expecting is the number of 0s and 255s.

The answer is posted as a reply to Otsu's method problem? . So it is now a duplicate I guess.


Intuitive understanding of frequency I know that $f = \frac{\text{cycles}}{\text{unit of time}}$  and I know that period $T = \text{time it takes to complete one cycle}$.

So if we take the definition

$f = \frac{1}{T} = \frac{\text{1 cycle}}{\text{time it takes to complete one cycle}}$

So if we have some $x$ cycles we get that $f = \frac{x \text{ cycles}}{\text{time it takes to complete one cycle}}$.

How can we have $x > 1$ if it takes $T$ to complete a cycle? In other words, how could $x > 1$ have completed in the time it takes to complete one cycle?

I’m guessing that $T$ is somehow decoupled from $x$ and I’ve got my definitions mixed up.

It's about units and dimensions. And I'm probably not the guy to cast an answer on that but I would like to state it as follows. In short; length, time, mass etc. are dimensions, whereas meter, seconds, kg are units of them respectively in SI. (for other dimensions and units, a physics book might help.)

Now, if a wheel rotates 7 times in 2 seconds; it's rotation frequency $f$ is found as:
$$ f = \frac{7 ~\text{ cycles}}{2 ~\text{ seconds}} = 3.5 ~\text{cycles-per-second (cps)}.$$

As can be seen, frequency is given by a (cycle) number divided by time. The cycle number has no physical units or dimensions; it's just a repetition count. Whereas the duration time has a unit of seconds (in SI).

Thus the dimesion of frequency is 1/time, whose SI unit is $1/s = s^{-1}$. This unit is also, and more frequently, denoted as Hz where it (implicitly) stands for (cycles) per second indeed. So, $3.5$ Hz stands for $3.5$ cycles per second.

Now, similarly the period is defined as the duration time per one cycle. Again note that since cycle is a number and does not posses a unit, then the dimension of the period is time and its unit is seconds in SI. But you can always denote it also as seconds per cycle, rarely practiced though.

Finally, using both notations you have
$$ T = \frac{1}{f} $$

So if frequency $f$ is defined with a unit of Hz ($s^{-1}$), then the period will have a unit of seconds (s) in SI. Otherwise if you define frequency as cycles per second then the period unit will be seconds per cycle.

So if you define the frequency as $3.5$ cycles-per-second , then the period associated is:
$$T = \frac{1}{f} = \frac{1}{ 3.5 \frac{\text{cyc}}{s} } = \frac{2}{7} \frac{\text{s}}{\text{cyc}}$$

Otherwise, in accepted norm, as was you would define freqquency as $3.5$ Hz. its period in SI will be:
$$ \boxed{ T = \frac{1}{f} = \frac{1}{ 3.5 s^{-1} } = \frac{2}{7} ~\text{ s} } $$


Question about SNR in continuous and discrete domain How to make sure the discrete domain's SNR equals to the corresponding SNR of continuous domain?

  If the definition of SNR in continuous domain is:

  $$SNR=\frac{\int_T|s(t)|^2dt}{noiseVariance} $$

  if we sampled $s(t)$ ($T_s$ is the sample period) then what is the
  corresponding SNR of the discrete domain? Is

  $$SNR=\frac{\sum_{n=1}^{N}|s(nT_s)|^2dt}{noiseVariance} $$

  ?

In your $snr$ expression, if you put a $1/N$ ratio in front of the numerator summation you'll have the standard expression for SNR with regard to sampled data (digital data).


How to Improve the Kalman Filter for Tracking the Periodic Motion of a Car? I have a quite typical Kalman filter to design. I really read a lot of articles about the design of this filter but the performances of my filter are still quite bad.
Here is my situation. I have a small car that does periodic constant motion. This little car is moving uniformly in a straight line with v in time t, with a velocity of 0 in time 4t (You can think of it as the little car moving at a constant speed v to a target place and do something, and then moving at a constant speed v to the next location.) I can measure the speed v of the car (with noise) and I use the model as follows to estimate my little car:
X_ =  X_last
P_ = P_last +Q
Here is my Matlab code: (I don't want to use the Matlab Kalman function ;) )

clear
clc
t = 1:628;
Z = square(t/50,50);
Z = Z+1.4;
Z = Z*80;
ZZ = randn(1,628);
ZZ = ZZ * 10;
Z = Z+ZZ;
plot(Z);

data = zeros(1,628);

Q = 0.1;
R = 3;
X_last = 0;
P_last = 1;
Dx = 0;

for i=1:628
    X_ = X_last;
    P_ = P_last +Q;
    Kg = P_/(P_+R);
    X = X_+Kg*(Z(i)-X_);
    P = (1-Kg)*P_;

    Dx = X_last;
    P_last = P;
    X_last = X;

    data(i) = X;
end
plot(t,data,t,Z);

My question is, is there a more appropriate model of a Kalman filter for the type of car that I'm trying to predict? Although the model works well for constant velocity, there's a trailing when the velocity goes from V to zero as Fig. shown. Is there a good solution to that? Thanks!
PS:
I control the velocity of the car by the current of its electrical motor, so I used the motor current to do matlab calculations（ y axis）

You basically have 4 models here:

Accelerating to constant speed.
Moving at constant speed.
Decelerating to zero speed.
Standing.

So the basic solution is building the 4 models and switching using Hard Switch between them.
Yet there is a smoother framework to handle smooth transition between them called Interacting Multiple Model (IMM) Kalman Filter.
Using the IMM Framework (Which is basically adding weighing step for the models) with the knowledge of the Input as described in other answers will give you a much better results.


Obtain the frequency response of a time domain chirp signal I have simulation data in the time domain and I need to represent them in the frequency domain.

The data have been obtained time integrating the following equation:

$m \ddot{x} + c\dot{x} + kx = A \sin{(\phi(t))}$

where $\sin{(\phi(t))}$ represents a logarithmic chirp. Therefore, since also the integrated time histories are chirps, I cannot use the FFT to transform them in the frequency-domain, but I have to do something like it was done in this answer: system response: time vs frequency. Why do I get different magnitudes?

Therefore I used the spectrogram in this way:

% Resample at a constant sampling frequency
dt = 0.0001;
t0 = time(1);
tend = time(end);
new_t = t0:dt:tend;
Fs = 1/dt;

acce_TS = timeseries(acce_cog,time);
acce_res = resample(acce_TS,new_t);

% Spectrogram
win = hamming(256);
noverlap = 0;
nfft = 512;
X = abs(spectrogram(acce_res.data, win, noverlap, nfft,Fs));
X = 2*X/sum(win);
surf(X)

And here come the problems:

1) On the time axis I have wrong values: the time history lasts $49.5 \ s$ whereas in the spectrogram plot I see values from 0 to 2000

2) The peak response should be at nearly $30 Hz$ (which is the resonance frequency of the system, given as $f_n = \sqrt{k/m}$), whereas from the spectrogarm I see it at nearly $5 Hz$

3) The amplitude of the peak of the spectrogram is a bit higher than in the time history

4) I actually need to see only the frequency-amplitude plot, not the 3D representation, so how can I do that? At the moment I am using view([90 0]), but is there a way to plot directly the 2D plot?

EDIT: Here the plots

1) On the time axis I have wrong values: the time history lasts 45 s whereas in the spectrogram plot I see values from 0 to 2000

RTFM on surf.  You'll need to do something like surf(new_t,new_f,X) instead of just surf(X).  Your new_f should be something like

new_f = linspace(0,Fs/2,256);

NOTE: I don't have access to matlab, and am not sure of the orientation of X. You may have to swap the parameter locations of new_t and new_f on surf.

  4) I actually need to see only the frequency-amplitude plot, not the 3D representation, so how can I do that? At the moment I am using view([90 0]), but is there a way to plot directly the 2D plot?

Do something like:

plot(new_t, max(X));

You may have to do max(X') instead, depending on the orientation of X.

I'll need to think about answers to 2 & 3. Watch this space.


Detect fixed overlay on set of images This is quite an "applied" question, but I don't see any more appropriate StackExchange, so here it goes. My situation is the following:

I have a sizeable (> 40) set of images (photos) that have been scanned using a flatbed scanner.
Unfortunately, there was a "smudge" due to some item previously placed on the scanner. This smudge remained fixed for all the scans, since it was due to some small amount of translucent residue on the glass bed.
The photos may be rotated some multiple of 90 degrees and offset, but are otherwise exactly what was scanned (high-resolution TIFF images).

Would it be possible to extract the "smudge" (i.e., the equivalent of an image filter applied to the physical photo that results in the scanned image), and then apply the inverse of this filter to all my images (with appropriate rotation and offset)? I don't know if tools already exist that could do this, but any guidance as to algorithms/approaches would be much appreciated.

Would it be possible to extract the "smudge" (i.e., the equivalent of an image filter applied to the physical photo that results in the scanned image), and then apply the inverse of this filter to all my images (with appropriate rotation and offset)?

It is difficult to reply Yes or No to this question. The results depend very much to what you are working with and for what purpose.

There are two things that are happening here:

The smudge absorbs some light
The smudge refracts and scatters some light

Dealing with 1 is easier than dealing with 2. Light is absorbed and therefore the smudge leaves behind it a shadow. If you had a "mask" of how much light passes through each point of the smudge, then you could get rid of it by adjusting the image pixel values (i.e. the comment suggestion).

If the smudge is caused by a deposit (e.g. finger markings, droplets, gel, etc) then it not only has diminished the light that went through it, it also has changed its path. Therefore, even if you correct the shadow, the "focus shadow" will still be there.

It is exactly like having an extra lens on top of the existing focusing mechanism.

This second bit could be corrected with a "filter", but it's easier said than done.

The main problem here is that you will not be able to recover an exact "image" of the smudge, so that you can then invert it. This is even worse if the smudge has been changing slightly with every image acquisition.

Short of obtaining an "empty" scan, you could try to obtain the average frame over N (grayscale) images. This is assuming that the N images have the smudge at exactly the same position. If not, then you will have to align them.

If N is high enough (at least 20 but it depends on image content too) then you will start to see the smudge shape over a (mostly) average background.

This can be used to diminish the shadow (at least) or deconvolve the affected regions to diminish the "focus shadow" too.

But, whether this is successful also depends on the use of the images. If you had to do measurements on the affected pixels, or look for patterns, this filtering operation could distort your results.

  I don't know if tools already exist that could do this, but any guidance as to algorithms/approaches would be much appreciated

These operations can be applied with low level image processing tools (e.g. Octave or skimage ) but for a quick test you could use:

Hugin to align the images
Gimp to stack the (aligned) images (in layers) and obtain the average image that would reveal the smudge.

The deconvolution bit is a bit more involved. There is ImageJ which has an astrophotography plugin that would be very similar to your use case except that in astrophotography the "smudge" is assumed symmetric and spatially invariant (same across the image)....which yours might not be. There are however other plugins you might find useful in this.

Hope this helps.

(If a subset of the images could be shared, we could have a quick go at it).


Image processing on a surveillance image I was given this question on a job interview and think I really messed up. I was wondering how others would go about it so I could learn from this experience.

  You have one image from a surveillance video located at an airport which includes line of people waiting for check-in. You have to assess if the line is big/crowded and therefore additional clerks are necessary. You can assume anything that may help your answer. What would you do?

I told them I would try to:

Segment the area containing people from the rest by edge detection
Use assumptions on body contour such as relative height/width to denoise unwanted edges
Use color knowledge; but then they asked how to do that and I didn't know

Thanks for sharing your experience. It's a bit off-topic, but my my first advice for a job interview would be to check the company's technology and products. This usually helps a lot in finding answers.

In this general case of crowd surveillance in a closed place, I would check the following items first:

get a reference image of the empty location in order to be able to detect people's location (even inaccurate) by background subtraction. This gives a rough estimate of the queue length;
check if the camera is facing the line, in which case I would run a face detector to get an approximate count of people in line. This provides a second estimate;
if necessary, learn (from the detected faces) the skin color in the image, then again segment the image based on skin color.

There are more advanced techniques available, but I think it would be overkill here.


Fastest implementation of fft in C++? I have a MATLAB program that uses fft and ifft a lot. Now I want to translate it to C++ for production. I used OpenCV but I noticed that OpenCV's implementation of fft is 5 times slower than MATLAB's. Then I tried armadillo but it was even slower. It was 10 times slower than MATLAB. Now I wonder is there any implementation of fft in C++ that is fast enough to compete with MATLAB?

Matlab's fft functions are all based on FFTW (this is confirmed here), so I guess the obvious choice for you should be FFTW. FFTW is hardware-independent but it can take advantage of some hardware-specific features.


map two RGB images in Matlab I wonder is there any algorithm or method that I can map two RGB images together in Matlab. I mean assume I have following images

Is there any way I find what number I should multiply each channel of one them to get the other one, or any other method.

I have tried something like: looking at destination image and see what color is dominating then adjust the other image by multiplying by different value using brute force method but that is very inefficient and tedious.

It is easier to do this in HSV color space than RGB color space.
Assume the first image(the pinky one) is imported into matlab workplace as variable A,and the second(the yellowish one) as variable B.
First, I removed the white margin of both A and B with a same Rect:

Rect=[82.51 32.51 511.98 509.98]; % I found this with [~,Rect]=imcrop(A);
A_c=imcrop(A,Rect);
B_c=imcrop(B,Rect);

I am going to find a simple way to transform image A_c to B_c in HSV color space;

A_c_hsv=rgb2hsv(A_c);
B_c_hsv=rgb2hsv(B_c);

Read the wikipedia page of HSV.. you will find that most different parts  between the two images must be in the Hue parts,i.e,  between A_c_hsv(:,:,1)  &  B_c_hsv(:,:,1) ;

And the average difference between their hues is:

ave_dif_in_hue=mean(mean(B_c_hsv(:,:,1)-A_c_hsv(:,:,1)));

I just used the hue parts in my transform:

A_c_hsv_transformed=A_c_hsv; % make a copy of A_c_hsv
A_c_hsv_transformed(:,:,1)=A_c_hsv_transformed(:,:,1)+ave_dif_in_hue; % change the hue parts only

let's check our transform:

figure;
subplot(1,3,1);imshow(A_c);title('original A_c');
subplot(1,3,2);imshow(B_c);title('original B_c');
subplot(1,3,3);imshow(hsv2rgb(A_c_hsv_transformed));title('transformed from A_c');

here is the result:

you can optimize the transform further, with dealing with the other two channel in HSV color space..
Or.Maybe something more needed.


Real Time Goertzel Algorithm Why is Goertzel Algorithm considered a block algorithm? Given that my input is bounded, couldn't I just run it forever (taking every sample that comes out after some length N) given a big enough word size for the intermediate coefficients?

I understand that I would get rectangular window spectral leakage, which may not be desirable, but let's assume it is good enough. If my input signal is stationary, then eventually I should get fairly good spectral resolution, correct?

But what if it isn't stationary? What is the memory of the Goertzel algorithm? How fast does the input signal frequency need to change before I get garbage results?

You can use forms of the Goertzel algorithm as a fixed length block filter, as a IIR filter with a very long impulse response, or as a fixed length running filter.

As a fixed length filter, a Goertzel results in a FIR filter similar (almost identical) to the result of 1 bin of a DFT or FFT of the same length.

As an IIR filter, a Goertzel filter is on the edge on instability (pole right on the unit circle), so you may need to bound the filter state and/or output externally, or reduce the gain slightly below 1.0.

As a fixed length running filter, you subtract the output of  2 offset (but stabilized) Goertzel IIR filters to get a FIR filter or 1 bin of a DFT (length equaling the offset) that updates every input sample, and with very few MACs per sample.


Group delay using cross correlation? An input signal + (thermal) noise travels through a channel and can be measured at its output.

Is it possible to use cross correlation between input and output signal to measure group delay? If not, is there a better way to do this?

You can use the cross correlation to measure the time delay between a signal and a delayed version of itself.

You can use cross correlation to estimate  the transfer function.

You probably would benefit the most from PHAT

https://www.mathworks.com/matlabcentral/fileexchange/26256-calculate-ambiguity-function-and-utilities?focused=5143816&tab=function


Can harmonics be created or simulated a la carte I am trying to learn about adding harmonics (audio) with a goal to make a plugin that can allow adjustment of the amplitude of each harmonic being added. I would not have the ability to write this but wanted to get an opinion about whether it is possible.  Example: to add a 2nd harmonic and adjust its level separately, add a 3rd and adjust its level etc. Ideally have a GUI that allows the ability to add harmonics and turn a knob or slider to adjust the amplitude of each harmonic up to around the 10th. The goal is to find out from some someone if it is possible and if so I’d pursue someone to write this for me.
Since the input can be any type of Audio from 20-20k, filtering is not likely the idea. I understand to some limited degree the methods of adding harmonic distortion, but have not been encouraged towards a la carte control. The end result is to custom tailor a harmonic distortion effect.
Thanks

Without re-hashing all of the maths from RBJs response, I made a crummy little plug in that did more/less what you are proposing. It had a stupid simple GUI with one slide for nominal amplitude and several more to control harmonic levels. It was OK, but not as useful as one would hope for audio work. The idea was that given a sinusoid with a nominal input amplitude, it would calculate the polynomial required for the given harmonics.  The bummer is that audio doesn’t work that way, at all.  The sound doesn’t have any more harmonic ‘richness’ if you can call it that, it just sounds distorted in different ways.
To wit, my peripheral understanding of products advertising harmonic distortion use crossover filters to divide the signal into different bands, on which distortion is applied. By narrowing the bandwidth, you get something more closely approximating a sinusoid. Again though, my impression was that they were fine, but not exactly what I had in mind.
My experience was that as simple as ‘adding harmonics’ seems, it’s a pretty tough thing to implement. My apologies for not providing a more rigorous analysis, but figured some first hand experience might be relevant.


Least Squares Re Synthesis of Short Time Fourier Transform (STFT) Working with my own MATLAB implementation of the short-time Fourier transform (STFT), I've managed to write code for the analysis step where a 1D time-domain signal $s[t]$ is progressively windowed, taken into the Fourier domain and arranged in a 2D matrix $S[t,\omega]$.  Each element in the 2D matrix is a function of time $t$ and frequency $\omega$.

This presentation is a rather nice overview of the STFT, and gives a number of equations detailing the analysis and re-synthesis steps.

However, I would like to be able to arbitrarily modify $S[t,\omega]$, and then use the re-synthesis to get back $s[t]$.

I believe that I should be able to change $S[t,\omega]$ in whatever way that I want, and then obtain $s[t]$ by the re-synthesis procedure.  This seems to be very similar to the idea of the phase vocoder.

As noted in section 3.1 of the presentation, the time-domain signal $s[t]$ can be recomposed using a least-squares procedure.  This is given as Equation (6) in the 1984 paper by Griffin and Lim.

The least-squares procedure is required to be applied when $S[t,\omega]$ is modified in some way.

Question:

What does Equation (6) of the Griffin and Lim paper mean?
What steps do I follow to numerically implement Equation (6)?

In the presentation, the equation is written in a slightly different way:

$$x(n)=\frac{\sum_{l=-\infty}^\infty w\left(n-lI\right)y\left(lI, n\right)}{\sum_{l=-\infty}^\infty w\left(n-lI\right)^2}$$

Note that $x(n)$ is the re-synthesized time-domain sequence, $w(n)$ is the window function, and $y(n)$ is the time-domain version of a column of the 2D matrix.

Steps:

From the presentation, here are the steps that I think are required to do the re-synthesis:

Let w_n be the discrete window vector and y_n(:,k) be the time domain vector computed using the IFFT on a column k of the 2D matrix.  Both w_n and y_n(:,k) are the same length.
Then, using Matlab syntax, we compute the point-by-point multiplication:
w_n .* y_n(:,k)

Is this the numerator of the expression above?
What happens during steps 3 and 4?
What do the infinite summations signify?

The least-squares re-synthesis procedure is very similar to the overlap-add (OLA) procedure.

Let w_n be the discrete window vector and y_n(:,k) be the time domain vector computed using the IFFT on a column k of the 2D matrix. Both w_n and y_n(:,k) are the same length.

Then, using Matlab syntax, we compute the point-by-point multiplication with the window:

w_n .* y_n(:,k)

As mentioned in the comments above, the y_n is trimmed to be the same length as the window w_t.  The w_n .* y_n(:,k) is then overlap-added in the same fashion as shown in the code associated with my previous post on the inverse STFT here.  The overlap-added sequence is the numerator of the expression shown in the original question above.

The same overlap-add operation is applied to the squared window w_t.^2.  This is the denominator of the expression shown in the original question above.

Then, the final re-synthesized output is simply the point-by-point division of the numerator by the denominator.  In Matlab syntax,

numerator ./ denominator


Planar Homographies in MATLAB I have taken an image in height 2 meters of the ground as the original image and I want to simulate another image in 50 cm's of the ground, such as the simulated image.

The focal length of my camera is 3.7 mm (f = 3.7).

I have used from planar homography: H = k'*(eye(3) - ((t*n')/d))*inv(k) for each plane in original image in MATLAB where in the camera coordinate system:

x axis is perpendicular to buildings in the right side of the ground,
y axis is parallel to the ground,
and hence z axis is upward perpendicular to the ground.

The origin is camera location in original image and k is 3x3 camera calibration matrix:

[f 0 0;
  0 f 0;
  0 0 1]

t is the transfer matrix [0 0 -1.5], n is the normal vector of a plane, d is the plane distance to the origin, and eye(3) is camera rotation.

For example, the plane of the ground has n=[0 0 1] and d = 2, is it true?

Please help me with implementation of other planes (sky, buildings) in MATLAB.

Original image:

Simulated image:

No can do. Your scene is not not planar, nor the camera motion you want to simulate is a pure rotation about the camera's center (a "pan-tilt"). These are the only cases where an image transformation with a pure homography will yield a realistic result.

One thing you can do is select one of the planes apparent in the image, and then
warp the whole image applying a homography to it - for example, you could make the left-side wall appear frontal. But your example has real parallax between the views, so it cannot be achieved with a homography.


Contourlet transform for image compression Before today I had no idea there is any transform an such; I was only aware with DCT and wavelet transforms but today one of my teacher who is doing research on image processing project had told me that contourlet transform is way better than wavelet transform.

It is better because wavelet transform extract frequencies in 3 directions (diagonal, horizontal, vertical) while contourlet transform is capable of extracting information in any direction.

There are various books on wavelet transform but I couldn't be able to find a book that mainly focus on contourlet transform, is there any book you would recommend for contourlet transform?

P.S: I'm not talking about research papers.

I do not know of any book solely devoted to contourlets. There is a chapter in Beyond Wavelets, 2003 (G. Welland, Ed.) by Do and Vetterli, named Contourlets. You also have a 100-page paper (almost a book) where contourlets are put in the perspective of Multidimensional Filter Banks and Multiscale Geometric Representations (2012).

At the risk of mildening your teacher's assertions, good compression is the result of a transformation with nice qualities and an efficient coding scheme. The "best" transform can fail to properly compress an image if you do not have a neat coding scheme associated with it. Contourlets, for instance, are slightly redundant. And its compaction capabilities are thus limited by the redundancy. To my knowledge, redundant transforms have not proved yet they are much better than critically-sampled ones for compression (except perhaps at large compression ratios).
Even if you are not talking about research papers, let me mention a couple of different contourlet versions. The CRISP-contourlet, which is non-redundant, and The Nonsubsampled Contourlet Transform, fully redundant.

When trying to understand a transformation, especially with wavelet-like decompositions, it is often useful to look at the interplay between the math and the sampling schemes.

In addition to Minh N. Do thesis, I am trying to find a copy of A. L. da Cunha PhD.


Time-invariant and Time-varying Systems Determine whether the following system is time-invariant or not:

$y(t)=x(t)\sin 10\pi t$

Solution:

Given: $y(t)=x(t)\sin 10\pi t$

$y(t)=T[x(t)]=x(t)\sin 10\pi t$

The output due to input delayed by $T$ sec is:

$y(t,T)=T[x(t-T)]=y(t)|_{x(t)=x(t-T)}=x(t-T)\sin 10\pi t$

The output delayed by T sec is:

$y(t-T)=y(t)|_{t=t-T}=x(t-T)\sin 10\pi (t-T)$

$y(t,T) \ne y(t-T)$

Conclusion: The delayed output is not equal to the output due to delayed input. Therefore, the system is time invariant.

  My query/doubt:

When the output delayed by T sec :

$y(t-T)=y(t)|_{t=t-T}=x(t-T)\sin 10\pi (t-T)
=x(t-T)\sin  (10\pi t-10\pi T)$

In case the value of $T$ is $\frac{2n}{10}$; $n\ge 0$, then the above equation becomes equal to $y(t,T)$

ie,

$y(t,T) = y(t-T)$ , When $T=\frac{2n}{10}$; $n\ge 0$

So can we say that the system is time-invariant when $T=\frac{2n}{10}$; $n\ge 0$ ?

Your conclusion is wrong: the system is time variant.

Regarding your question: for a system to be time invariant, it needs to be time invariant for any delay.

Even if there are specific delays when the system behaves as time-invariant, in general the system still is classified as time-variant.


How does amplitude scale when FFT size is smaller then signal duration? I have one quite simple question: how does the amplitude of the signal scale when we do FFT where FFT window length is shorter then duration of the signal? Does this means that if FFT size is 12 time shorter that also amplitude of the signal will scale by 1/12? if so what is the mathematical explanation.
Thx in adavance for any reply.

When you perform simple DFT of your signal with length $N$, then you have to scale it properly. Normally it is done by dividing magnitude by length of spectrum vector, which is $N$. Another thing is that normally we are interested in one-sided spectrum, therefore you need to multiply your above result by 2 to compensate energy. This will assure, that values of magnitude are representative.

In case when you are using windows different from rectangular (Hamming, Blackman, etc.) you are not dividing by $N$ - you must scale it by sum of your window samples.

Therefore when you are windowing your signal, you get a new signal with a length of a window. Thus in your case amplitude will be scaled down by 12.

Example for 1Hz sine wave and two window types is shown below. Please note, that frequency axis is chopped at 3Hz. You can see that amplitude of peaks is equal to unit amplitude of signal. Also in case of Hamming window side-lobes are much lower than in case of rectangular, sacrificing main-lobe width at the same time. This can be seen better in logarithmic scale.

MATLAB code used to produce result:

clc, clear all, close all
f0 = 1;        % fundamental frequency
fs = 50;       % sampling frequency
T = 6*pi;      % upper time limit
t = 0:1/fs:T;  % time vector
N = length(t); % number of samples

% 5 Hz sinusoid with unit amplitude
s = sin(2*pi*f0*t);
% Windowing signal
win = hamming(N)';
s_win = s.*win;
% Magnitude spectra
S = abs(fft(s));
S_win = abs(fft(s_win));
% Frequency vector
f = (0:fs/N:round(fs/2));
% Take one-sided spectrum
S_h = 2*S(1:length(f));
S_win_h = 2*S_win(1:length(f));
% Normalisation
S_h = S_h / length(s);
S_win_h = S_win_h / sum(win);
% Time domain signal plots
subplot(2,1,1)
plot(t,s)
xlim([0 T]); grid on; hold on;
plot(t,s_win,'r');
legend({'Rectangular window','Hamming window'})
xlabel('Time'); ylabel('Amplitude')
% Frequency domain plots
subplot(2,1,2)
plot(f, S_h)
grid on; hold on;
plot(f, S_win_h, 'r');
legend({'Rectangular window','Hamming window'})
xlabel('Frequency'); ylabel('Amplitude')
xlim([0 3])


Is the sum of white noise and shifted white noise white noise again? Let $W[k]$ be a stationary white noise with variance = 1

Question: Is $X[k] = W[k] + c \cdot W[k-1]$ white noise?

$c$ is a real number.

Calculate the autocorrelation of the process.

$$\begin{align}
R_{xx}[n] &=\mathbb{E}[(W[k] + c W[k-1])(W[k-n] + c W[k-1-n])] \\
&=\mathbb{E}[W[k]W[k-n]]+ \mathbb{E}[cW[k]W[k-1-n]]+\mathbb{E}[cW[k-1]W[k-n]]+\mathbb{E}[c^2W[k-1]W[k-1-n]] \\
&=\sigma^2\delta[n]+c\sigma^2\delta[n+1]+c\sigma^2\delta[n-1]+c^2\sigma^2\delta[n]\\
&=\sigma^2(1+c^2)\delta[n]+c\sigma^2\delta[n+1]+c\sigma^2\delta[n-1]
\end{align}$$

The definition of white noise implies that $R_{xx}[n]=\sigma^2\delta[n]$, which is not the case here.


Is there a widely available implementation of an adaptive, recursive, numerically stable IIR filter? I need to identify the coefficients of a linear, causal, time-invariant physical system that can be described by a classical state-space formulation.

For the sake of the example, suppose that the system has two states, only one of which is observed, and one input. Suppose that this state-space is described by the following matrices:

  a =
               x1          x2
   x1   -1.25e-05    1.25e-05
   x2   1.389e-05  -2.083e-05

  b =
              u1
   x1          0
   x2  6.944e-06

  c =
       x1  x2
   y1   1   0

  d =
       u1
   y1   0

The discrete, equivalent, zeroth-order-hold transfer function is given by:

  3.481e-05 z + 3.446e-05
  -----------------------
   z^2 - 1.97 z + 0.9704

In filter terminology, this is clearly an IIR filter because of the higher-order coefficients in the denominator.

I'm looking for a suitable adaptive filter algorithm that will:

be numerically stable
allow for online updates, i.e. be recursive

From my research so far, I've read that if you treat the recorded signal values as inputs to this filter, then it's possible to treat this filter as a FIR filter. This is called the equation error approach. One can then apply, for example, the QR-RLS adaptive FIR filter to find the coefficients.

But in my trials on the real system I often find filter coefficients that are not physically possible, and I just don't know if my hand-crafted adaptive filter is correct or not. I would like to try with an "official" implementation of such an adaptive filter.

Therefore, my question is whether there exists (for MATLAB or in C) a freely available implementation of an adaptive IIR filter that has the same advantages of a QR-RLS FIR adaptive filter.

While I may not be the final authority, I have looked into this. In short, no.

From my own research, the problem is stability. When the filter is extended to have a wider bandwidth, which of course means that the y[n] values and lags 'move quickly', and then is quenched to a lower bandwidth, where the y[n] values and lags 'move slowly', the y[n] buffer may be loaded with fast oscillated entries that appear as initial conditions to the low bandwidth filter. The filter can easily go unstable. If the bandwidth is altered in an adiabatic manner this instability may be avoided (I'm supposing) but then there's an additional constraint on your adaptivity.

Now, you don't have to take my word for it. In "Adaptive Filters: Theory and Applications" (2nd) page 5 (Amazon), Boroujeny writes, "However, as we shall see in the later chapters, because of the many difficulties involved in adaptation of IIR filters, their application in the area of adaptive filters is rather limited." The author states that the adaptation process can place (digital) poles outside of the unit circle, even if they started within the unit circle. The rest of his book is mostly dedicated to FIR filters.

I hope this information helps.


Perform lowpass filter by hand I've done research and could not find any resources on this topic.

I have message signal composed of more sinusoids. Then I preform amplitude modulation, that is multiply that signal with higher frequency carrier signal, and obtain modulated signal. Then to demodulate that signal, I multiply the signal with same carrier signal, and obtain demodulated version of signal. Then I should pass the signal through low pass filter, to obtain signal identical to message signal.

If I do this by matlab using fdesing.lowpass it is fine I get good output, but how to do this low pass filter by hand on paper, if I just ignore frequencies higher than the maximum frequency in message signal I don't get same amplitudes.

I multiply signals by using inverse Euler formula for $\cos$.

Any link or help would be appreciated.

If I didn't make myself clear, let me now and I will try to write an example (which I omit for now because I don't see here any latex support for writing mathematical equations).

how to do this low pass filter by hand on paper?

You can use very simple low pass filter - Moving average. Filter is quite simple - it is sum of N sample of signal divided by N. If length of filter will be equal to   period of carrier signal then performance of this filter will be good.


Reason of Fast Fourier and Inverse Fast Fourier I am not very good at mathematics
I was doing some image processing so I came to know about FFT and IFFT
I was learning about them but the mathematics involved seems to be too dense
So can someone tell me about what are FFT and IFFT and their general application in short in LayMan's Language
So that again when I go to learn it I get it a little easier.
thanks

FFT and IFFT are algorithms that implement the (inverse) discrete Fourier transform (DFT). These transforms convert a signal into another representation, namely the frequency domain, and back. This other representation allows us to analyze certain properties of the signal and also to change these properites which would be hard to do in the original representation.

Consider a digital audio signal: it's a sequence of samples that represent amplitude over time. This representation is useful to study the temporal properties of the signal and to play it back to make it audible. We can hear that the audio signal contains low and high sounds (frequencies) but from looking at the signal it's impossible to tell what frequencies exactly. After applying the DFT to the signal we can analyse what frequencies it contains and how strong they are. We can even amplifiy or attenuate certain frequencies and then apply the inverse DFT (IDFT) to make the signal audible again. In this way we have built a filter. The application of IDFT is necessary in this case, because the frequency domain representation of the audio signal is useless for playback.

Image processing deals with two-dimensional signals and therefore two-dimensional (I)DFT. The concept of frequencies with images is not so illustrative as with audio signals. Edges, i. e. sudden changes of intensity are represented by high frequencies whereas smooth gradients are represented by low frequencies. This allows us to do edge detection, for example. If high frequencies are attenuated this will result in a smoothed image after application of IDFT. Again, for just looking at the image, the frequency domain representation is useless but it allows analysis and filtering.


Spectrum analysis : how to represent FFT results I'm trying to program a spectrum analyser in C++ but I'm faced with a conceptual problem : how to map the results of the FFT to an image ?

Let me detail :
I'm doing all this in a dynamic way.

I have chunks of audio data as float [-1;1] arriving to my class.
I put them in a circular buffer which has the same size as my FFT window
I apply a Hann(ing) windowing function
I apply an FFT to the result.

This part works great.

Now the result of my FFT is a two dimensional array (my library automatically does the sqrt(re^2 + im^2) part) with bin indexes and amplitudes.

My goal is to draw this result from time to time into a window. The problem is, the window has way less pixels than there are bins (if my buffer size is 65536, I have 32768 bins), so I have to do some kind of "averaging".
I could do a mean, a root square mean, take the maximum, etc etc ... I have no idea what is the best method to use ? The application of this will be music, more specifically, detect frequency masking by plotting several semi-transparent sepectrums on the top of each other.

Any advice on a good way to do that would be much appreciated

Actually, what do the squared magnitude FFT result represent? They represent the total energy in a given frequency bin. So to downscale your results to less bands, you can sum them. It's not even necessary to average the results - however it would have the same effect, as you map the results to an image, thus rescaling them.

However, note that if you need to rescale by something else than an power of two, you'll need to do it using supersampling, as the results won't map directly to your new bins.


discrete sine transform using FFT i have made a subroutine to calculate 1d fft of a series of size 2^n. I need to calculate discrete sine transform of a series using the fft subroutine. I am not sure what to do.

well, you can do a little Googling and Wikipeding and find the definitions of the DST and pick one of them.  The get the definition of the DFT (the FFT is an efficient method to compute the DFT), use Euler's formula on the DST definition and see where a DFT pops out.

i might get back to this tonight and spell it out, but if you're not sure where to go next, the above (hit the links) tells you where to approach the problem.


What is the maximum data rate a wireless link can support using BPSK for given BER and SNR? What is maximum data rate a wireless link can support using BPSK for 10^-6 BER, Bandwidth of channel is 200 kHz and SNR 21.335dB ?
As per my understanding , Eb/No can be calculated for BPSK for given BER using curves and value found is 11.29dB, and after that we can use relation

Eb/No = S/N * W/R

R comes out to be around 2 Mbps! Is it correct? Need your input!

In BPSK, BER is equal to $Q(\sqrt{2E_b/N_0})$, so you need $E_b/N_0>11.3$ to have a BER no larger than $10^{-6}$. Since the bandwidth is 200 kHz, the maximum rate is 400 kb/s. So, since $SNR\approx136$,

$$\frac{Eb}{N_0}=\frac{S}{N}\frac{W}{R}=136\frac{200\cdot10^3}{400\cdot10^3}=\frac{136}{2}=68$$

Since $68>11.3$, you can be sure that you can transmit at 400 kb/s and meet your BER requirement.

Note that I'm assuming baseband transmission (so that the bit rate is twice the bandwidth). Note also that the question, as posed, does not make a lot of sense, since the bandwidth will essentially determine your bit rate.


Which come first, PCA, FFT, filtering Might be a simple question, but I keep second guessing my algorithm, so a confirmation would be good. I first do bandpass filtering, then PCA, then FFT of the top principle components.

I also do some sine wave fitting on the filtered data, and on the principle components, but it doesnt make sense to do it on the fft data? Right?

Thank you for your patience on this stupid question.

Functionally PCA, FFT could be interchanged without changing anything.

e.g.

[U,~,~] = svd(A);
[FU,~,~] = svd(fft(A));

fft(U) and FU are the same (up to an arbitrary phase rotation)

However, it may help speed to do BPF first, then FFT, and discard the stopband to reduce dimensionality before doing PCA.


ROC and impulse response For the LTI system given below, there are three regions of convergence.

$$H(z)=\frac{5-3z^{-1}}{1-\frac53z^{-1}-\frac23z^{-2}}$$

a) Find all possible regions of convergence for this filter.

b) For each region of convergence, find the impulse response of the filter.

I tried this problem using the poles and came up with the following outcomes. Want to confirm

Are the ROCS $z<2$ and $z>-\frac{1}{3}$  and $z=3/5$? Any help in this regard is highly appreciated.

It seems that you've correctly found the two poles: $z_{\infty,0}=2$ and $z_{\infty,1}=-\frac13$. The only other thing you need to know is that the possible ROCs are limited by the pole radii, and that there can be no poles inside an ROC.

Note that $z<2$ doesn't make much sense because $z$ is a complex number. If you mean $|z|<2$ then the pole at $z=-\frac13$ would be inside that ROC, so $|z|<2$ is not a valid ROC for the given transfer function.

The three possible ROCs for your example have the following forms:

$|z|<a$
$a<|z|<b$
$|z|>b$

The numbers $a$ and $b$ are of course related to the poles of the transfer function. I trust that you can take it from here.


Calculating bpm from machine sound (not music) I have a hydraulic breaker on a excavator that hammers with a dynamic bpm. I have recorded this sound so that it is in a mp3-file and would like to determine the momentary bpm for the length of the file.

I have tried using processing with beat detection from the minim-library but so far without any good results.

Is using beat detection with energy mode in a selected band the preferred way to go? I have had troubles getting good information about where, in terms of frequency, the beats are.

An example of the sound could come from http://www.youtube.com/watch?v=596Dtm95Txo

There are two differences between your problem and music beat detection:

The beat in music is a "latent" grid. It's not because a track is at 120BPM that you have exactly 2 notes per second. In music, the beat is more of a lowest common periodicity than an exact pulse train - that's why beat trackers often make octave errors, and also why there is no consensus between human listeners as to what is the BPM of some tracks. Your task doesn't suffer from the same ambiguity because every "hit" counts.
The beat in music rarely changes over time; and when it changes, it is often progressive. As a result, the tempo can be assumed constant over windows of 5 or 6s for music. Judging from the video, you cannot make this assumption; because the "rhythm" here is very irregular; with pauses/breaks, and frequent changes.

Because of these, I fear that all the assumptions embedded into beat tracking algorithms made them useless on your problem. My suggestion would be to just use an onset detector with a high threshold (to make it more robust to noise). The events you want to detect seem to cover a rather wide frequency band anyway and have a fast decaying envelope.

Then, simply measure the time delta between adjacent events (when this delta is too high; consider that this is a pause); and eventually, average this over groups of a couple events.


Why is there a negative in front of the phase response equation for this complex exponential? first time on here!

I'm working through "Digital Signal Processing using MATLAB" by Vinay and Proakis. Good book.

I am stuck on this example tho.

Shouldn't the imaginary part in the denominator (magnitude calculation) be (-0.9sin(w))^2 since the Frequency Response equation has an e^-jw?

Also why does the angle calculation include a -arctan instead of just a regular arctan.

Thanks

For the magnitude, the denominator term will be $$1-0.9e^{-j\omega}=1-0.9(cos\omega - jsin\omega) = 1-0.9cos\omega+0.9jsin\omega$$

Hence the plus sign.

For the phase, if you multiply the numerator and denominator with the complex conjugate,

$$H(e^{j\omega}) = \frac {1-0.9cos\omega-0.9jsin\omega}{(1-0.9cos\omega+0.9jsin\omega)(1-0.9cos\omega-0.9jsin\omega)}$$

The denominator will become real.
Therefore,
$$\angle H(e^{j\omega}) = arctan(\frac {-0.9 sin\omega}{1-0.9cos\omega})$$
Taking the minus sign outside is how you get the negative in the phase.


Phaser effect with feedback control IIR filter calculation This question is possibly relevant to this one, but it's not duplicate as far as I'm aware.

I'm trying to make phaser with feedback control as shown in the picture below: (Got it from here)

First problem for me was that feedback loop here doesn't have any delay, so at time $n$ I should pass allpass cascade output $y_n$ to its input at this exact time which seemed impossible to me. However, thanks to this comment I was able to figure out how to transform IIR filter cascade to incorporate this feedback into it in terms of its coefficients:

Recursive equation for $k$-th order IIR filter:
$$y_n = b_0x_n + b_1x_{n-1} + ... + b_kx_{n-k} - a_1y_{n-1} - a_2y_{n-2} - ... - a_ky_{n-k}$$
As we pass system output through explicit feedback loop our input sequence $x_n$ is defined as follows:
$$x_n = d_n + fy_n$$ where $d_n$ is system "dry" input and $f$ is feedback coefficient such as $|f| < 1$

If we plug this into the first equation we'll get
$$y_n = b_0(d_n + fy_n) + b_1(d_{n-1} + fy_{n-1}) + ... + b_k(d_{n-k} + fy_{n-k}) - a_1y_{n-1} - a_2y_{n-2} - ... - a_ky_{n-k}=\\
=b_0d_n + b_0fy_n+...+b_kd_{n-k}+b_kfy_{n-k}-a_1y_{n-1}-...-a_ky_{n-k}=\\
= b_0fy_n+b_0d_n+...+b_kd_{n-k}-(a_1-b_1f)y_{n-1}-...-(a_k-b_kf)y_{n-k}$$

Moving all $y_n$ terms to the left and grouping them we'll finally get our desired system coefficients:
$$(1-b_0f)y_n = b_0d_n+b_1d_{n-1}+...+b_kd_{n-k}-(a_1-b_1f)y_{n-1}-(a_2-b_2)y_{n-2} -...-(a_k-b_k)y_{n-k}$$

Now, the problem is: we all know we cannot directly use filters with order $k > 2$ without stumbling upon various stability problems. Say, I want to implement phaser with four notches, meaning I'll have allpass cascade of overall order $k = 8$ In that case I would probably want to find poles/zeros of my system analytically (instead of coefficients) and split it to second order sections afterwards or something like that, but I cannot figure out how to find poles/zeros instead of coefficients or is it even possible to do so?

All I've got is that using the fact that for allpass of order $k$ and having poles $[p_1, p_2, ..., p_k]$ its zeros will be always equal to $[1/p_1, 1/p_2, ..., 1/p_k]$ and knowing that all poles of initial system (without feedback) are equal I can write down my problem as follows:

$(x - p)^8 + a(x-1/p)^8= 0$

$p,a \in ℝ$

$|p| < 1$

$|a| < 1$

Also, below you can find matlab code which calculates filter coefficients for the whole phaser, you can see unstable case on
[bf, af] = notch_reso2(100, 48000, 0.5, 0.5);

function [ b, a ] = notch_reso2(  fc, fs, feedback, depth  )
order = 8;
w = pi * fc / fs;
a1 = (1 - cot(w)) / (1 + cot(w));
k = a1 ^ order;
p = -a1;
z = 1/p;

zs = repmat(z, 1, order);
ps = repmat(p, 1, order);

b = poly(zs) .* k;
a = poly(ps);

f = feedback;

af = a - b .* f;
bf = b ./ af(1);
af = af ./ af(1);

bf = bf .* depth + af * (1 - depth);
b = bf;
a = af;
end

Solution can be found in the well-known paper of Vadim Zavalishin. In sections (4.1) and (4.2) there's an explanation on how to implement zero feedback loop around cascade of LP filters. He even mentions later in the paper that exactly the same approach can be used to implement phaser (6.1)


ML estimation - solve for x I'm trying to solve the following maximum likelihood estimation but for multiplicative noise instead of additive noise:
So the goal is to do ML-estimation for a scalar constant $x$, which is multiplied with noise:
$$\tilde{d}[n] = a[n]x$$
The noise $a[n]$ is a realization of the i.i.d. random variable $A \sim \mathcal{N}(\mu_A, \sigma_A)$.
The probability density function is given by
$$ \pi(x,\mathbf{\tilde{d}}) = \frac{1}{x^N} \pi_A\left(\frac{\tilde{\mathbf{d}}}{x}\right)$$
with $\pi_A$ being the pdf of $A$.
Now I took the natural logarithm of this expression and derived with respect to $x$, which let to the following term:
$$\frac{\partial}{\partial x} \ln\left(\pi(x,\mathbf{\tilde{d}})\right) = -\frac{N}{x} + \frac{1}{\sigma_A^2} \sum_{n=0}^{N-1} \left[ \frac{\tilde{d}[n](\tilde{d}[n]-\mu_a x)}{x^3} \right] \overset{!}{=}0$$
But how do I continue from that? My idea was to solve this equation for $x$, but I can't figure out a closed form solution...

It looks okay to me. If you define empirical mean $\hat{\mu}_d = \frac 1N \sum \tilde{d}[n]$ and empirical second moment $\hat{\gamma}_d = \frac 1N \sum \tilde{d}^2[n]$ then you effectively have an equation of the form $$-N/x + A/x^3 + B/x^2 \stackrel{!}{=} 0,$$ where $A$ and $B$ directly depend on $\hat{\mu}_d, \hat{\gamma}_d, \mu_A, \sigma_A$. If you multiply this equation with $x^3$ it becomes a quadratic equation in $x$, which you can solve for $x$...
I get something like $$x^2 - \frac{\hat{\mu}_d \mu_A}{\sigma_A^2} x + \frac{\gamma_d}{\sigma_A^2} = 0$$ but that was scribbled down in a rush so I would have to check it again later. Can you confirm or correct?


Adaptive LMS Algorithm MATLAB I'm having some trouble implementing my LMS Adaptive Filter in MATLAB to separate wideband and narrowband signals from a voice signal.

I'm using a delayed version of my input as a reference as well as the error term.

step = 0.01;
w = zeros(1, N);
xDelayed = [zeros(1, 100) x']'; % delaying input

for n=1:length(x)
    e = x(n) - w(1)*xDelayed(n);
    w = w - step*e*xDelayed(n);
end

It's essentially an implementation of this

$$w(n+1) = w(n) - \alpha e(n) x(n)$$

For some reason, my entire w (N long) vector is all the same value.
UPDATE:

M = 5;
N = length(sound)
w = zeros(M, N);
STEP_SIZE = 0.01;
d = sound;
x = sound_delayed(1:N);

for i=(M+1):N
   e(i) = d(i) -  x((i-(M)+1):i)*w(:,i);
   w(:,i+1) = w(:,i) + mu * e(i) * x((i-(M)+1):i)';
end
for i=(M+1):N
    yd(i) = x((i-(M)+1):i)*w(:,i);
end

There are several problems in your code. First, it looks like you're confusing iteration and vector indices. The computation of e should use all values of the current (delayed) data vector, filtered with the current filter coefficients. In the update equation, you subtract a scalar from a vector, which is not what you want. Again you should be using all values of the current data vector (the length of which must equal the chosen filter length).

Take a look at the Matlab code in this question (the second one in the EDIT-part). I haven't run it but it looks like it deals correctly with the vectors in the update loop.


Chroma key algorithm with unknown background color Suppose that I have an image with a foreground subject and a monochrome background. The image is a real world picture so the pixels belonging to the background do not share the same exact value. I want to isolate the foreground but I don't know the background color in advance. What are the state of the art algorithms to achieve this goal?

Something self-written should be easy. "Background" probably means that the majority of pixels at the border of your image have that color.

So, just get all the border pixels, put them in an array, convert to hue and then do a histogram. Use the maximum histogram bin as your chroma key color.

Starting with border pixels of that color, erode your image on a hue gradient, i.e. make sure that the hue doesn't change too much from its neighboring pixels, and then add pixels to your background mask until you're done OR
just do classical chroma-keying and convert all pixels of colors "close enough" to that color to background mask.


Compensation for Irregular Time Step for DFT (FFT) I want to run a clustering algorithm (svm, knn) on the ferquency spectrum data of a temperature sensor that published at irregular times.
Here is the temperature data to take the FFT:

I got the mean frequency for this temperature data as well as every temp sensor in the dataframe the average is every 18mins, longest distance between 2 points is 3hrs:

# group by short id
groups = df.groupby('id')
time_diff = groups.apply(lambda df: df.published_at.diff().mean())

Isolated the graphed temp data in a series:

signal = df.loc[df['_id'] == 'A1']['temperature']

Stored size of signal and mean sampling frequency as variables:

# sampling frequency:
Fs = time_diff[:1]
Fs
Out[217]: 00:18:54.085526

# size
S = signal.size
Then took the fft and calculated dBs

X = np.fft.fft(signal)
X_db = 20*np.log10(2*np.abs(X)/S)
And plotted the results:

plt.plot(X)
plt.show()

plt.plot(X_db)
plt.show()

These graphs intuitively do not look like they correspond to the original data. The objective with the fft is to then classify the data using SVM, however, I am not sure which variant fft is appropriate, nor if using the mean time frequency is either.

One option, in addition to the Non-uniform DFT that has been mentioned, is simply interpolating the data on a regularly-spaced grid by using a timestep approximately the same size as the smallest timestep, and then taking the DFT on the interpolated data.


gaussian noise in an image according to
Estimating acquisition noise P. 2

Estimating the mean and the standard deviation for each pixel is calculated  as below:

$$
     \begin{equation*}
      i, j = 0,\dots,N-1 \\
      \bar{I}(i, j) = \frac{1}{N} \sum_{k=0}^{N-1}  I_k(i,j)\\
      \sigma(i, j) = \left [ \frac{1}{N-1} \sum_{k=0}^{N-1}(I_k(i, j) - \bar{I}(i,j)^2) \right]^{1/2}
     \end{equation*}
$$

My Question

What is $I_k$ ?
What should mean of each pixel be? Is it the mean of whole Image?

Thanks

To estimate quantities about a pixel at location $(i,j)$, you need a set of pixels in a neighborhood, close enough to $(i,j)$ in behavior, to compute statistics on.

Assuming that neither the camera nor the object at pixel $(i,j)$ move, the index $k$ may denote a sequence of images in time, as answered by @MarcusMüller.

From the other slides at CMPE 264:   Image Analysis and Computer Vision, it is not evident that the presenter considers image sequences. A  possibility is that the index $k$ in $[0,N-1]$ denotes a subset of $N$ pixels somehow around $(i,j)$, maybe both in time or space.

In space, $I_k$  could denote the  $k$th pixel  from a subblock of the whole image with $N$ pixels total, centered around the pixel with coordinates $(i,j)$. Such a lousy notation could avoid cumbersome notations for the subblock that moves with position $(i,j)$.

Pixels around the pixel with coordinates $(i,j)$ are considered as realizations of a random process modeling the center pixel, and serves to provide an estimate of the average of  this pixel.


Pattern classifier experiment for determining sex from face image of persons using kNN? I am very much interested in pattern classification experiments and its theory etc. So I am trying some home brew experiments using MATLAB. As a first major experiment, I would like to classify persons based on sex from facial images. My doubt is basically is it worth attempting it, because right now I don't have idea on what are the basic facial features to extract and how to do that? I want basic facial features to distinguish between males and females. I am using kNN classifier.

At your first stage, which classifier to use is not that important, you may need to make your code run first. There was a very famous paper on Face Recognition and Gender Determination, in which a very complex transform was implemented to represent each node on the face. The used gabor kernal based wavelet transform to get set of coefficients for each node(refers to jet), and use sparse jets to represent the face. They also consider the template matching to ensure that each face can be decomposed of the jets. They built a look-up table to record the correlation between each jet and the gender. When a new image comes, the gender is determined by the pooling result of the jet sets.

But most approaches nowadays won't be that complicated. Basically the image intensity is used as the face input (or RGB color components but they are intensity equivalently). Each image matrix is reshaped into a long vector, and all the images from training set will form a large M X N matrix where M is the image pixel size and N is the sample number. Then Principal Component Analysis is used on the matrix to extract the core components, and project them onto a sub-space. Intra-class and inter-class scatter difference is calculated, then Linear Discriminant Analysis is applied to maximize the ratio of the determinant of Inter-class scatter and intra-class scatter. Please refer to this link to see more details (The algorithm is implemented by Python though).


DFT symmetry vs DFT duality in Richard Lyons' "Understanding DSP" I am reading Richard lyons, understanding dsp, chap 3.

Article 3.2 is about property of dft symmetry but any where in this chapter, i am unable to find discussion about dft duality property

I want to know here ,author means both dft symmetry and dft duality to be a same thing?

1. Duality between discrete frequency and discrete time domain.

DFT Duality is generally referred to the duality of DFT-IDFT pairs. This in turn comes from the similarity between analysis and synthesis expressions of DFT and IDFT.
$$X[k] = \sum_{n=0}^{N-1}x[n]e^{-j\frac{2\pi}{N}nk}$$
$$x[n] = \frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{j\frac{2\pi}{N}nk}$$
(Although I will encourage you to understand DFT as projection of Finite-length sequence x[n] onto Discrete Fourier Orthogonal Basis Vectors, and IDFT as representation of x[n] as linear combination of those orthogonal fourier basis vectors, where coefficients of linear combination is are the DFT coefficients $X[k]$. You might want to read this answer : DFT as Projection on Orthogonal Fourier Basis Vectors)

By Duality between the DFT-IDFT pairs I mean, for example, DFT of discrete $\delta[n]$ is a constant $1$ in freq-domain and DFT of constant $1$ will be discrete $\delta[k]$ in freq-domain. And other example would be DFT of discrete $rect$ sequence i.e. sequence of $M$ $1$'s in an $N>M$ length sequence. It has its DFT as follows:
$$X[k] = e^{-j\frac{\pi}{N}(M-1)k}\frac{\sin[\frac{\pi}{N}Mk]}{\sin[\frac{\pi}{N}k]}, \forall \ k \in \ \{0,1,2,...,(N-1)\}$$
And, if we took DFT of $x[n] = e^{-j\frac{\pi}{N}(L-1)n}\frac{\sin[\frac{\pi}{N}Ln]}{\sin[\frac{\pi}{N}n]}$, we will get a sequence of $L$ consecutive $1$'s in an $N$ length DFT. We can just run the steps backwards to compute the DFT of this $x[n]$.

For continuous time fourier transform, duality would mean the following : if $f(t)$ has fourier transform $F(\Omega)$ then a time domain function $F(t)$ will have its fourier transform as $2\pi f(-\Omega)$.

  2. Symmetry in Discrete Fourier Transform Coefficients.

DFT Symmetry in the book is mentioned in context of symmetry properties in DFT coefficients when DFT of a real valued time domain sequence $x[n]$ is computed. That is DFT coefficients of all real valued $x[n]$ are conjugate symmetric modulo $N$.
$$X[k] = X^*[(N-k) \ \mod \ N]$$


Increase power of signal in Matlab I have a scenario where I have two OQPSK modulated signals and I need to increase the power of one signal by 3 dBm while keeping the other signal at the same power level. Could someone guide me on how to achieve this in Matlab.

To increase the power of x by 3 dB:

x .*= 10^(3/20)


How do the motion vectors work in predictive coding for MPEG? In MPEG, there is a process where an image is broken up into macroblocks and a motion vector is computed for each of those macro blocks.  You then transmit these vectors, along with the prediction errors, to reconstruct the next image in the video sequence.

I'm trying to get a firm grasp on how this works.  Each macroblock has a motion vector associated with it, which (if vector is [1,0]), says all the pixels in this block move 1 in the x direction and 0 in the y direction for the next frame.  If all the motion vectors don't align correctly, won't this leave areas of the image unaccounted for (like the area where that macroblock was in the first place)?

For example, I have the following question I found.

Consider the following image at time t:

7   7   7   7
7   7   5   5
7   5   5   8
8   8   8   8
9   9   9   9
9   9   9   9

This image was broken up into 2x2 macroblocks and the following motion vectors were sent along to recreate it:

(0,0)  (0,0)  (0,1)  (-1,1)  (0,0)  (0,0)

The image at the previous time step, t - 1, looked like this:

7   7   7   7
7   7   5   4
7   7   7   7
7   5   8   8
8   9   8   9
9   9   9   9

What were the errors transmitted?

How would you solve this?

To simplify your confusion - there two processes :

1. Motion estimation
2. Motion compensation

Before we talk about the estimation, we should talk about the Motion compensation.

Let say, the $ Image_{t}(x,y) $ is split in block $ Blocks_{t}[k](x',y') $.

The task of Motion compensation is to produce $ Blocks_{t}[k](x',y') $ from any region of
$ Image_{t-1}(x,y) $.

Hence another block not necessarily aligned at 16x16 boundary is a best possible match $ Blocks_{t-1}[k](x'+mx,y'+my) $

Here, $ mx, my $ is called motion vectors.

We can calculate the error between the target and reference as

$$ Err_{t}[k](x,y) = Blocks_{t}[k](x',y') - Blocks_{t-1}[k](x'+mx,y'+my) $$

So now, the encoder basically transmits $Err_{t}[k](x,y)$ (with DCT and quantization) and $ {( mx, my) }[k] $  for each block,.

So encoder has 2 work to do:

1. Motion Estimation
The process or estimating $ { mx, my }[k] $ for every $ k $ such that $ Err_{t}[k](x,y) $ is minimized is called Motion estimation.

2. Generation of error image after Motion Compensation
The process of constructing $ Blocks_{t}[k](x',y') $ from $I_{t}$ image pixels and $ { (mx, my) }[k] $ is called Motion compensation.
The error image is what gets transmitted.

Finally, decoder can re-do motion compensation on its' own using moiton vectors and the error image to make final re-construction of the image.

Now we realize a few points:

Best Motion estimation helps minimize the energy that is required to be transmitted and hence optimizes the bits for a given quality.
However, even if $ { (mx, my) }[k] $ is not ideal or if the scene has significant change over last picture, the $ Err_{t}[k](x,y) $ is always transmitted to receiver - hence the reconstruction is always perfect (modulo the loss created by the quantization). Hence, even if you have suboptimal motion vector or redundancy is not much, the reconstruction is still always perfect albeit with more bits!
Every block $ Blocks_{t}[k](x',y') $ is motion compensated on its own right - hence even if the actual motion vectors of any neighboring blocks has no effect in the construction. Hence, it is not necessary to have motion vectors perfectly aligned for perfect reconstruction to be possible.
Albeit algorithm  exists who are smart enough to guess that if $Blocks_{t}[k]$ has motion vector $( mx ,my )[k]$ the guess for $Blocks_{t}[k+1]$ might just be closer to that only.
Finally suppose the next picture is completely different it is possible that energy for $$ Energy (Err_{t}[k](x,y))  > Energy ( Blocks_{t}[k](x',y') ) $$.

In such cases, it might be more advisable to transmit block directly without prediction than to send the difference. This is also possible in the encoder by a provision called INTRA block.


Which mother wavelet for a scalogram? I'm attempting to create a real-time scalogram (from a 1 dimensional signal) in the style of a spectrogram;

Looking through various papers+books; the Gabor wavelet, or complex Morlet seems to be favoured for keeping a close relation to frequency.

Though I was hoping to use a real valued wavelet though due to computational complexity concerns... What wavelet would be recommended?

The mother wavelet of your scalogram should have a similar shape to the usual peak shapes you want to detect (I suppose you use it to detect peaks of your signal).
However, I would like to ask you what would you like to use wavelets for? I could give you a more specific answer for your question.


Complexity of FFT derivation I am confused regarding the complexity of the Fast Fourier Transform (FFT). The Discrete Fourier Transform is:

$$\qquad X\left [ k \right ]=\sum_{n=0}^{N-1}x[n]W_{N}^{kn}\quad \text{where}\quad W_{N}=e^{\frac{-j2\pi}{N}}\tag{1}$$

The first step in deriving the FFT is to express $(1)$ as the sum of even and odd inputs

$$\qquad X[k]=\sum_{n\ \rm even}x[n]W_{N}^{nk}\ +\ \sum_{n\ \rm odd}x[n]W_{N}^{nk}\tag{2}$$

And after many other derivations the final equation $(3)$ is obtained as

$$\qquad X[k]=\sum_{r=0}^{N/2-1}x[2r]W_{N/2}^{rk}\ +W_{N}^{k}\ \sum_{r=0}^{N/2-1}x[2r+1]W_{N/2}^{rk}\tag{3}$$

Equation $(3)$ consists of $2$ length $N/2$ FFTs, which is more computationally efficient than performing one big length $N$ FFT.

My question is: isn't equation $(2)$ also $2$ length $N/2$ FFTs? In other words, why all the bother about deriving equation $(3)$? Isn't equation $(2)$ already more computationally efficient than equation $(1)$? Or perhaps I'm missing something?

Any help is greatly appreciated.

No, Equation (2) is not the sum of 2 FFTs (of lengths $N/2$) as you claim, and it is not computationally more efficient. The canonical way of evaluating Equation (1) is via Horner's method:

$$X[k] = (((\cdots (x[N-1]\cdot\omega + x[N-2])\cdot\omega + x[N-3])\cdot\omega + \ldots)\cdot\omega + x[0]\tag{a}$$

where $\omega = W_N^k$ which requires $N-1$ multiplications and $N$ additions: lather, rinse, and repeat for for other values of $k$. Equation (2) just splits this Horner's rule computation into two parts with no change in complexity. That is, we compute $\sum_{n\ even}x[n]W_{N}^{nk}$  as $$(((\cdots (x[N-1]\cdot\Omega + x[N-3])\cdot\Omega + x[N-5])\cdot\Omega + \ldots)\cdot\Omega + x[0]\tag{b}$$ where $\Omega = \omega^2$ and
$\sum_{n\ odd}x[n]W_{N}^{nk}$ as
$$(((\cdots (x[N-2]\cdot\Omega + x[N-4])\cdot\Omega + x[N-6])\cdot\Omega + \ldots)\cdot\Omega + x[1])\cdot\omega.\tag{c}$$ Each computation requires $N/2$ multiplications and $N/2$ additions with no savings in computation over the standard brute-force Horner's method applied to a vector of length $N$ in $(a)$.

Equation (3), however, is the core of the FFT butterfly; the two sums in Equation (3) are DFTs of two shorter length sequences.  Define \begin{align}x_{\text{even}}[k] &= x[2k], &k = 0, 1, \ldots, N/2-1,\\
x_{\text{odd}}[k] &= x[2k+1], &k = 0, 1, \ldots, N/2-1,
\end{align} and let the DFTs of these shorter sequences be denoted by $X_{\text{even}}[k]$ and $X_{\text{odd}}[k]$ respectively. Then the computation of Equation (1) has been re-written in Equation (3) as
$$X[k] = X_{\text{even}}[k] + W_N^kX_{\text{odd}}[k] \tag{**}$$
which requires one multiplication and one addition per bin. Note that $k$ varies from $0$ to $N-1$ in $(**)$, and the shorter DFTs in $(**)$ are assumed to be periodic functions of period $N/2$ (where's rb-j when I need him?)  But, now note that $X_{\text{even}}[k]$ and $X_{\text{odd}}[k]$ are DFTs of length $N/2$, and so, if $N/2$ is an even number (hint: this means that $N$ is a multiple of $4$), we can use the idea behind $(**)$ to express $X_{\text{even}}[k]$ as
$$X_{\text{even}}[k] =  X_{\text{even, even}}[k] + W_{N/2}^kX_{\text{even, odd}}[k]$$
and $X_{\text{odd}}[k]$ as
$$X_{\text{odd}}[k] =  X_{\text{odd, even}}[k] + W_{N/2}^kX_{\text{odd, odd}}[k]$$
that is, both DFTs of length $N/2$ have been expressed in terms of 2 DFTs of length $N/4$, and so on.  If $N = 2^m$, we can proceed in this way to where we have $N/2$ DFTs of length $2$ that we compute using one multiplication and addition per DFT bin. So, $N$ multiplications and additions total to compute these $N/2$ DFTs. (Purists whining that a length $2$ DFT needs only one addition and one subtraction are kindly requested to go jump in the lake!). Next, we combine these $N/2$ DFTs of length $2$ into $N/4$ DFTs of length $4$ using one multiplication and one addition per DFT bin for a total of $N$ multiplications and additions to do the combination. Next, we combine these $N/4$ DFTs of length $4$ into $N/8$ DFTs of length $8$ using one multiplication and one addition per DFT bin for a total of $N$ multiplications and additions to do the combination..... Ultimately, we arrive at $(**)$.  This iterative procedure for computing the DFT of length $N$ is (one of the many versions of) the FFT algorithm for computing the Discrete Fourier transform. It has $m=\log_2N$ stages requiring $N$ multiplications and $N$ additions at each stage for a total of $Nm = N\log_2N$ multiplications and additions. In contrast, the use of Horner's rule would need $N(N-1)$ multiplications and $N^2$ additions. Hence the soubriquet FAST for this kind of algorithm.

Purists who have survived their jumps into the lake can now start their cacophony about radix-2 versus radix-4 versus radix-8 transforms and the exact number of multiplications and additions needed for implementing an FFT algorithm.


Efficient real-time pitch detection algorithm What is the most efficient pitch detection algorithms in terms of computation? I need something to implement on a DSP.
The typical signal in my situation would be a mixture of low-level wide-band noise plus a series of strong harmonics below 1000 Hz (think of the noise produced by industrial machinery or by car engines). The system needs to sample ("know") the current fundamental frequency of the signal every second or so.

YINS fft is the best method that worked for me in the past. I have my open source implementation of it written in awesome D language: link . I compared the results with Praat they match.


What does it mean when one model fit the data better than does a comptitor model? I'm trying to review MLE (maximum likelihood estimation).

What does it mean when one model fit the data better than does a comptitor model?

You need a criteria for saying whats better (i.e. a risk function). You say estimator A is better than estimator B with respect to a given risk function if the risk of estimator A is better than estimator B.

In the case of estimation, the risk is often the mean square error, i.e. if $\hat{X}$ is your estimator and $X$ is the thing you're trying to estimate, $R(\hat{X}) = E[\lVert X - \hat{X} \rVert_2^2]$ where $\lVert \cdot \rVert_2$ is the 2-norm. Then, you say estimator A is better than estimator B if it has lower mean square error.

When working in the mean square error framework for deterministic parameters, you have the notion of efficient estimators, i.e. ones which meet the Cramer-Rao lower bound, which are the best you can do in some sense (as described on the wiki article).


Adding phase noise to a signal results in amplitude variation on a constellation - whats the mistake? *plots inside* I am attempting to add phase noise to a modulated signal.

I am doing this by taking a phase noise mask (single sided 0 Hz to 1
MHz), which shows dBc values away from central value or reference
value. Phase noise is decreasing as we move in frequency away from
the reference point which is 0 Hz. [Top left image in the 4-plot
image below]
I then make this double sided.
I convolve this phase noise plot with a spectrum of a clean complex
exponential at 1 MHz, this is my oscillator for converting frequency.
I think this has imparted the phase noise onto the complex
exponential.  [Top right image in the 4-plot image below]

The spectrums are convolved  to make a oscillator with phase noise:
$\mathcal{F}(e^{j \omega t}) * \mathcal{F}( e^{j \phi(t) }) = \mathcal{F}(e^{j \omega t})* \mathcal{F}(e^{j B}) = \mathcal{F}(e^{j( \omega t + \phi(t))}) $

where $\mathcal{F}$ is fourier transform, $ \omega$ = 1 MHz, $\mathcal{F}( e^{j \phi(t) }) $ is the double sided phase noise mask on 0 Hz, then also this can be interpreted as a bandwidth $B$ of the upper and lower sidebands of the phase noise centered on 0 Hz (the double sided spectrum in the top left plot below when not in logarithmic x-axis scale). The convolution result is the top right image, symmetric phase noise translated from 0 Hz onto my a oscillator or CW or sine wave at 1 MHz.

I take the Inverse Fourier Transform of this convolution to produce
an oscillator in the time domain, my thought is that the phase noise
is on it.
I then multiply my complex modulated signal (16-APSK) by this oscillator (complex exponential with phase noise) to produce an up conversion and then downconvert with a clean complex exponential
with no phase noise. This should leave the phase noise on the
modulated signal.

The math now being done in time domain:
Up-conversion using the inverse transform of the result above
$A(m) e^{ j( \phi (m))}  e^{j( \omega t + \phi(t))} $
Down-conversion of the upconversion
$ A(m) e^{ j( \phi (m))}  e^{j( \omega t + \phi(t))} e^{-i \omega t}  $
where $A(m) e^{ j \phi (m)} $is a modulated signal

However I am showing results that are not like classical phase noise,
which would be rotation acrcs around constellation points but it
appears as an amplitude modulation?  [Final constellation image]

I know this is wrong, but I stuck why. If I subtract in the time
domain the clean sinusoid without phase noise and the oscillator with
phase noise, I get this variation across the full time domain
simulation. [Central images are the real and imag values of the oscillators in time and the bottom images are the difference]

Why is this producing amplitude variation and not phase variation?
Is there another way to do this?

--Edit Updated constellation based on replies below--

As the OP has described, the noise as shown is AM noise only as the constellation points are not changing in phase at all. A double sided spectrum with symmetric sidebands in phase with the carrier is amplitude noise:
Consider a reference carrier as 1 at angle 0 degrees for which we do the following to create sinusoidal sideband noise by adding to it a small sinuoid with with $\alpha$:
$$c_{n(am)} = 1 + \alpha cos(\omega_nt)$$
This is a real waveform which only changes in amplitude and all samples will be of zero phase (AM noise only). We can gain further insight by decomposing the cosine using Euler's Formula (of which the magnitude and starting phase of each are the two sidebands in the Fourier spectrum):
$$c_{n(am)} = 1 + \alpha/2 (e^{j\omega_nt} +e^{-j\omega_nt})  $$
However if the two sidebands were in quadrature relationship to the carrier, then for small $\alpha$ where the small angle approximation holds, the noise will be dominantly phase noise:
$$c_{n(pm)} = 1 + j\alpha cos(\omega_nt)$$
$$= 1 + j\alpha/2 (e^{j\omega_nt} +e^{-j\omega_nt})$$
For the case of noise density where the noise is relatively small so that the small angle approximation holds, the OP can simply multiply the sidebands representing the noise pattern by $j$, mapping the result to be PM instead of AM. Note from the graphics below that negating half the spectrum (for example multiplying all components in the negative half spectrum by $e^{j\pi}$ or 180°) will also accomplish converting AM to PM (to do this with an FFT use an odd number of bins, multiply all the $(N-1)/2$ upper bins by $-1$ and multiply the DC bin by j to have this come out perfectly):

For adding phase noise from a power spectral density in MATLAB or Octave, please see this great function with all that worked out from Alex Bur-Guy:
https://www.mathworks.com/matlabcentral/fileexchange/8844-phase-noise


LowPass Filter & Window *Img & Code* - Frequency Response has no lobe roll off Not sure if I am doing something silly but I am not seeing lobe roll off on the frequency response of my windows and the low pass filter I made in MATLAB.
Please can you help ... several hours just looking at this and my textbook -.-

It should be like this showing lobe roll off

%% Low Pass Filter
Fs=10*10^3; cutoff = 1000; %100 Hz
fc = cutoff / Fs ;
b = 800/Fs ;%Transition bandwidth of 800 Hz
M = 51 ; %Length of Filter
n = 0:1:M-1;
blackman_window = 0.42 - 0.50 * cos((2*pi*n)/(M-1)) + 0.08*cos( (4*pi*n) / (M-1));
hamming_window  = 0.54 - 0.46 * cos((2*pi.*n)/(M-1));
hanning_window  = 0.50 - 0.50 * cos((2*pi.*n)/(M-1));
blackman_window_ft = fftshift(fft(blackman_window)/M);
hamming_window_ft  = fftshift(fft(hamming_window)/M);
hanning_window_ft  = fftshift(fft(hanning_window)/M);

sinc_func = sin(2*pi*fc.*(n-((M-1)/2))) ./ (2*pi*fc.* (n-((M-1)/2)));
sinc_func((M+1)/2) = 1 ;
sinc_func_ft = fftshift(fft(sinc_func)/M);

h_lpf = sinc_func .* blackman_window;
h_lpf = h_lpf ./ sum(h_lpf); %Normalise
H_lpf = fft(h_lpf);
H_lpf = fftshift(H_lpf);

figure(1);cla;clf
subplot(3,2,1);hold on;
plot(blackman_window,'DisplayName','Blackman')
plot(hamming_window,'DisplayName','Hamming')
plot(hanning_window,'DisplayName','Hanning')
plot(sinc_func,'DisplayName','Sinc')
legend
title('Windows and Sinc Impulse Response');hold off

subplot(3,2,2);hold on;
plot(20*log10(abs(blackman_window_ft)),'DisplayName','Blackman')
plot(20*log10(abs(hamming_window_ft)),'DisplayName','Hamming')
plot(20*log10(abs(hanning_window_ft)),'DisplayName','Hanning')
plot(20*log10(abs(sinc_func_ft )),'DisplayName','Sinc')
legend
title('Window and Sinc Frequency Response (dB)');hold off;

subplot(3,2,3);hold on;
plot(h_lpf);
title('Low Pass Filter Impulse Response');hold off;

subplot(3,2,5);hold on;
plot(abs(H_lpf),'DisplayName','Magnitude');
plot(real(H_lpf),'DisplayName','Real');
plot(imag(H_lpf),'DisplayName','Imag');
legend
title('Low Pass Filter Frequency Response');hold off;

subplot(3,2,6);
plot(20*log10(abs(H_lpf)));
title('Low Pass Filter Frequency Response (dB)');

You need to zero pad the FFT so that it reveals more of the frequency response (increasing the resolution of the plot). I corrected your code (while taking fft, I took 1024 point FFT for windows as well as h_lpf)

%% Low Pass Filter
Fs=10*10^3; cutoff = 1000; %100 Hz
fc = cutoff / Fs ;
b = 800/Fs ;%Transition bandwidth of 800 Hz
M = 51 ; %Length of Filter
n = 0:1:M-1;
blackman_window = 0.42 - 0.50 * cos((2*pi*n)/(M-1)) + 0.08*cos( (4*pi*n) / (M-1));
hamming_window  = 0.54 - 0.46 * cos((2*pi.*n)/(M-1));
hanning_window  = 0.50 - 0.50 * cos((2*pi.*n)/(M-1));
blackman_window_ft = fftshift(fft(blackman_window,1024)/M);
hamming_window_ft  = fftshift(fft(hamming_window,1024)/M);
hanning_window_ft  = fftshift(fft(hanning_window,1024)/M);

sinc_func = sin(2*pi*fc.*(n-((M-1)/2))) ./ (2*pi*fc.* (n-((M-1)/2)));
sinc_func((M+1)/2) = 1 ;
sinc_func_ft = fftshift(fft(sinc_func,1024)/M);

h_lpf = sinc_func .* blackman_window;
h_lpf = h_lpf ./ sum(h_lpf); %Normalise
H_lpf = fft(h_lpf,1024);
H_lpf = fftshift(H_lpf);

figure(1);cla;clf
subplot(3,2,1);hold on;
plot(blackman_window,'DisplayName','Blackman')
plot(hamming_window,'DisplayName','Hamming')
plot(hanning_window,'DisplayName','Hanning')
plot(sinc_func,'DisplayName','Sinc')
legend
title('Windows and Sinc Impulse Response');hold off

subplot(3,2,2);hold on;
plot(20*log10(abs(blackman_window_ft)),'DisplayName','Blackman')
plot(20*log10(abs(hamming_window_ft)),'DisplayName','Hamming')
plot(20*log10(abs(hanning_window_ft)),'DisplayName','Hanning')
plot(20*log10(abs(sinc_func_ft )),'DisplayName','Sinc')
legend
title('Window and Sinc Frequency Response (dB)');hold off;

subplot(3,2,3);hold on;
plot(h_lpf);
title('Low Pass Filter Impulse Response');hold off;

subplot(3,2,5);hold on;
plot(abs(H_lpf),'DisplayName','Magnitude');
plot(real(H_lpf),'DisplayName','Real');
plot(imag(H_lpf),'DisplayName','Imag');
legend
title('Low Pass Filter Frequency Response');hold off;

subplot(3,2,6);
plot(20*log10(abs(H_lpf)));
title('Low Pass Filter Frequency Response (dB)');