How do I calculate the range of a fixed-point number with $a$ integer bits and $b $ fractional bits?

Let's assume that we are dealing with unsigned number types. If you would use all $a+b$ bits for the integer part then the set of possible numbers would be:
$$\left\{0, 1, 2, 3, \dots, 2^{a+b}-1\right\}.$$

These numbers can be divided by $2^b$ (or multiplied by $2^{-b}$) to take use of $b$ bits for the fractional part, resulting in this set of possible numbers:

$$\left\{0, 1\times2^{-b}, 2\times2^{-b}, 3\times2^{-b},\dots, \left(2^{a+b}-1\right)\times2^{-b}\right\}
\\=\left\{0, 2^{-b}, 2\times2^{-b}, 3\times2^{-b},\dots, \underline{\underline{2^a-2^{-b}}}\right\}.$$

So your formula gives the largest number that can be represented (double underlined).

The formula can also be understood as going a step $2^{-b}$ or one least significant bit (LSB) worth backwards from $2^a$ which is the first number that has the same truncated binary string representation as 0, the first number in the system. In a similar way in the 8-bit unsigned integer system we take a LSB-sized step backwards from 256 (1 0000 0000 binary) to obtain 255 (1111 1111 binary) which is the largest representable number in that system.

I just found out that range has an established meaning in arithmetic: Largest value minus smallest value. The smallest value happens to be zero so in this case the largest value equals the range.


Why do I have to use variance and not mean value to calculate noise?

Typically the first moment (mean) is zero and therefore contains (relatively) no information. The most obvious thing then, to distinguish different levels of noise, is to look at the variance.

And if the variance happens to be zero, then surely there is no noise at all.


Discrete Kalman filter for a continuous system

The mass-spring-damper combination is an LTI system described by the following continuous-time linear differential equation

$$ \ddot{x} = - \frac{k}{m} x - \frac{b}{m} \dot{x} + \frac{1}{m} u $$

where $u$ is the deterministic input force (N),  $k$ is the spring constant (N/m) and $b$ is the damping coefficient (N.s/m).

Assuming two states as the position $x(t)$ and the velocity $\dot{x}(t)$; $X = [x , \dot{x}]^T$, and assuming a single measurement as the position $x(t)$ and assuming no process noise, and scalar deterministic input $u$ only added to the process equation, then the continuous-time linear dynamic system process and measurement equations for the Kalman filter can be written as

$$ \dot{X} = F X + G u $$
$$ y = H X + v $$

where
$$
F = \begin{bmatrix}
0 & 1 \\
-\frac{k}{m} & -\frac{b}{m}
\end{bmatrix}
$$

is the system dynamics matrix,

$$
G = \begin{bmatrix}
0 \\
1/m
\end{bmatrix}
$$
is the deterministic input matrix,

$$
H = \begin{bmatrix}
1  & 0
\end{bmatrix}
$$
is the measurement matrix, $v(t)$ is the measurement noise.

The first step in the derivation of the discrete Kalman filter is the discretization of the associated matrices by first finding the fundamental matrix $\Phi(t)$ which is given by :

$$ \Phi(t) = \mathcal{L}^{-1} \{ (s I - F)^{-1} \} $$

where $s$ is the complex variable of the Laplace transform.

For the given $F$, the fundamental matrix is found to be:

$$
\Phi(t) = e^{\frac{-b}{2m}t}
\begin{bmatrix}
\cos(\omega t) + \frac{b}{2m\omega } \sin(\omega t)  & \sin(\omega t) / \omega  \\
-\frac{k}{m\omega } \sin(\omega t) & \cos(\omega t) - \frac{b}{2m\omega } \sin(\omega t)
\end{bmatrix}\\
$$

where $\omega = \frac{\sqrt{4mk-b^2}}{2m}$.

Then the discrete fundamental matrix of the Kalman filter is found from $\Phi_k = \Phi(T)$ where $T$ is the sampling period.

$$
\Phi_k = e^{\frac{-b}{2m}T}
\begin{bmatrix}
\cos(\omega T) + \frac{b}{2m\omega } \sin(\omega T)  & \sin(\omega T) / \omega  \\
-\frac{k}{m\omega } \sin(\omega T) & \cos(\omega T) - \frac{b}{2m\omega } \sin(\omega T)
\end{bmatrix}\\
$$

We should also find the discrete versions of $H$ and $G$ which are found as $H_k = H = [1 , 0]$ and

$$G_k = \int_0^T \Phi(t) G dt .$$

By the help of a symbolic integration package (such as provided by MATLAB) one can compute the discrete deterministic input matrix $G_k$ as

$$
G_k = \frac{e^{aT}}{m (a^2+w^2)}
\begin{bmatrix}
\frac{a}{\omega} \sin(\omega T)-\cos(\omega T) + e^{-aT}  \\
 a (1-\omega) (\cos(\omega T)-1) + (a^2+ \omega) \sin(\omega T)
\end{bmatrix}\\
$$

where  $ a = -b/2m  $ and $\omega = \sqrt{4km-b^2}/2m$ as before.

After having found all the necessary discrete matrices $\Phi_k$, $G_k$, $H_k$ we can setup the discrete Kalman process and measurement equations through which the recursive solution can be applied.

$$X[n] = \Phi_k X[n-1] + G_k u[n]$$
$$ y[n] = H_k X[n] + v[n]$$

where $u[n]$,$v[n]$ and $y[n]$ are sampled versions of deterministic input, measurement noise and measurement respectively, and $X[n]$ is the sampled version of the state vector.

Then you will apply the Kalman recursion to estimate the states. Below are some simulation results from a MATLAB implementation of the associated Kalman filter. In the simulation, the true signal is generated by solving the associated ODE with a simple Runge-Kutta numerical method. Eventhough it's not required for this LTI system whose analytical solution is already available, it'll be helpful to use a numerical method later if more complicated (nonlinear) differential equations are encountered.

In this simulation, the model was exact. And therefore it's a little trivial. For more realistic applications, one should either replace the linear differential equation with more realistic nonlinear ones hence change the model and/or apply process noise to compansate for the lack of exact process model information.


Why must outputs be the same for a pair of inputs for a causal system?

So $x_1(t)$ and $x_2(t)$ are equal up to a certain time $t_0$, and diverge after that time. What this definition is saying is that, up to time $t_0$, the system will produce the exact same output, without regard for whether the input is actually $x_1(t)$ or $x_2(t)$.

If the system is non-causal, then it may "know" whether the input is actually $x_1(t)$ or $x_2(t)$, because it can "look ahead" in time beyond $t_0$, and it could modify its output accordingly, even before $t_0$.

As an example, consider a lightbulb connected to a battery and to a switch. The switch is initially on, and the lightbulb is also on. At time $t=0$, you decide whether you'll turn off the switch at time $t=60$ (one minute). The lightbulb, being causal, cannot know in advance what will happen, and will stay on at least up to $t=60$. A non-causal lightbulb will know in advance what will happen and may turn itself off at, say, $t=30<60$.

I hope this makes sense -- let me know if you need more clarification.


There is a problem in checking whether the homography is OK.

The algorithm for checking correct homographies may interest someone, so I will write it down here:

1) Create a quadrilateral $ABDC$ with vertex coordinates (in homogenous coordinates):

$$\begin{eqnarray} A:& (-w/2,-h/2, 1.0) \\ B:& (w/2,-h/2, 1.0) \\ C:& (-w/2,h/2, 1.0) \\ D: &(w/2,h/2, 1.0) \end{eqnarray}$$

where $w,h$ are width and height of the image, respectively. If the whole image frame (a rectangle) is transformed to convex quadrilateral, than any convex quadrilateral within it will also be transformed "neatly".

2) Create transformed quadrilateral $A'B'D'C'$ in which every vertex is transformed using the computed homgraphy (e.g. $C'=HC$ ). From now on, all points will be converted to non-homogenous coordinates.

3) Compute vectors $\vec{u}$, $\vec{v}$ for parametric representation of diagonals:

$$\begin{eqnarray}d_{1} :& A+(D-A)s =A+ \vec{u}s \\ d_{2} :& B + (C-B)t=B+\vec{v}t \end{eqnarray}$$

The intersection of diagonal comes from $d_{1}=d_{2}$:

$$t=\frac{1}{d}\left[(B_{y}-A_{y})\vec{u}_{x} - (B_{x}-A_{x})\vec{u}_{y}\right]$$

$$s=\frac{1}{d}\left[(A_{x}-B_{x})\vec{v}_{y} - (A_{y}-B_{y})\vec{v}_{x}\right]$$

Then the convex quadrilateral satisfies $s,t \in (0,1)$.

In practice, one can introduce a fudge factor to avoid not only non-convex and degenerate quadrilaterals, but also near-degenerate ones, like when three points are near colinear. So the test can be modified such that $s,t \in (\lambda, 1.0-\lambda)$, where lambda is a small number (e.g. $\lambda=0.01$).

Older problem, fixed in the above algorith:

I found the problem here - having a certain homography, the test can pass for a smaller quadrilateral, but not for the larger one. This is why some "ill" homographies passed through.

The green square represents a source image, the orange is a transformed one. As you can see, the left hand one is convex, but starts deforming as the source is larger:

Finally even larger source yield non conver quadrilateral:

I made a mistake with scaling. The points in homogenous coordinates $(x,y,w)$ were scaled in $x$ and $y$ direction, but in $w$. This is why the same transform produced different quadrilaterals.

I have corrected the algorithm accordingly.


Recommendations for learning and getting hands-on Real Time Signal processing

I had taken the course at my university. So most of the materials were provided as part of it. I worked with 2 different DSPs both of which were provided at the lab.
If you are looking to buy hardware, Real-Time Digital Signal Processing by Sen M Kuo and Bob H Lee is a good book with implementation and architecture descriptions of TMS320C55x. It has around 9 interesting application-level implementations as well. But getting all the hardware necessary might be too much though.


Proof from wave equation of upper limit of active noise cancellation

I don't have a mathematical proof for you, but here are 2 issues to think about.

The impulse response from the actuator (speaker) to the error measurement (microphone) is non-minimum phase. This is unavoidable because it takes a non-zero amount of time for an acoustic wave to travel from the speaker to the microphone. In a digital system there is additional latency introduced because you cant instantaneously sample, process, and output a signal. This latency limits your cancellation bandwidth.

The other issue is that even if you had a super low latency system and could get cancellation up to 10kHz it would still only be cancelling noise at the error microphone. At 10kHz your wavelength is about 3.4cm. The distance from the error microphone to the eardrum is a sizable portion of this wavelength, which means that whatever cancellation you thought you had at 10kHz is gone, and you might even be boosting now. So your perceived cancellation is limited by how accurately your microphone measurement approximates the sound field at the eardrum.


how to detect two parallel lines in a binary image with MATLAB

You mentioned the Hough transform, but your code doesn't use it. However, it can help you to find out the orientations of lines.

The maximum values of the Hough transform correspond to probable lines. These maxima are defined by their coordinates $(\rho, \theta)$, where $\rho$ is the distance from the origin to the closest point on the straight line, and $\theta$ is the angle between the $x$ axis and the line connecting the origin with that closest point (see Wikipedia).

Therefore, two parallel lines have a close $\theta$ (or $\theta + \pi$).


In the first case, $x(0)$ and $x(1)$ and $x(2)$ etc are all uncorrelated. (And since they're gaussian, they're independent). This means that, to get a realization of the first random process, you have to take $\infty$ independent samples from a gaussian distribution, and assign each of these $\infty$ samples to a $x(n)$.

The second process is much simpler. You just take one sample from a gaussian distribution (and call it $\alpha$), and then assign all of the $x(n)$ to this $\alpha$. In this way, $\alpha$ may be random, but you will get $x(0)=x(1)=x(2)=\cdots=\alpha$.

I'll post some plots, as you asked. I used the following code in MATLAB:

N = 3;
T = 10;
t = 0:T;
x1 = randn(N,T+1);
x2 = repmat(randn(N,1), 1, T+1);

figure
subplot(3,1,1); stem(t, x1(1,:));
subplot(3,1,2); stem(t, x1(2,:));
subplot(3,1,3); stem(t, x1(3,:));

figure
subplot(3,1,1); stem(t, x2(1,:));
subplot(3,1,2); stem(t, x2(2,:));
subplot(3,1,3); stem(t, x2(3,:));

In this first figure, I plot three independent realizations of Random Process 1:

And in this figure, three independent realizations of Random Process 2. Note that in this case, all the $x$ signals are constant, but the actual level changes from one realization to another:


Determining the 'amount' of speaking in a video

The most straightforward and easy thing you could do is rectify the signal (i.e. take the absolute value of its samples) and then run a moving average filter to get an approximate estimate of the "envelope" of your signal. Imagine the moving average filter like estimating the mean of the signal values between samples 0-100, then 1-101 then 2-102 and so on.

The result of this will be an inaudible slow varying waveform. You can then do a histogram on its values and use that to estimate a threshold which could then be used to "tag" the interesting parts of your video.

However, I feel that this is a lot of hard work for this kind of problem.

Therefore, I would also like to suggest that you use a video editor that presents the video "signal" as two separate ribbons, one for the actual video and another representing the amplitude of the sound. You can then quickly inspect the audio ribbon and estimate quickly the parts that could be seen as "interesting" and use the video editor to either directly transcribe those areas or cut the video down into smaller segments. Unfortunately, at the moment, I can`t recall a freely available video editor with this feature but if I do I will amend this response.


Difference Between Correlation / Convolution and Matrix Multiplication

Well I will try to explain.

Let us first discission in time domain:

1) Let us say you have two signals, x and y. By Convolution in time domain, you mean that you flip(invert) one of the signals(lets say ,y) and then slide it over the other(x). By corerlation, you just don't flip them, rest is the same.

2) Now when we talk about correlation, we first discuss 1-D case, where you have just one signal(x) and then for multiple signals(x,y...)

One Variable:

The variance for x will tell you how spread or varied that signal is, i.e. a small will correspond to signal which is changing very much.

As as engineers, we like that just at looking at value we can make prediction about the nature of the signal. For that reason, variance is not a good estimate, as you cannot compare signals easily by just looking at their variance. For that reason we introduce Standard Deviation, which is infact the underroot of variance, but by looking at it you can easily compare the nature of two signals.

Multi Variable:

As variance is in one signal case, in mutli variable, where you have X,Y ...etc signals, we replace it by Covariance. A covariance will tell you how much X and Y vary with each other. But again like variance, we cannot compare just based on co-variance, so we introduce correlation(Standard deviation in one signal case). Now by looking at correlation, you can compare the signals. And then a correlation coefficient,which is just the normalization of correlation, and it ranges from 1 to -1. If the correlation coefficient between 2 signals(x and y) in this case, near to 1, means these both signals are same or same behavior, whenever one signal goes up, the other also and vice versa. 0 means that they have no relationship at all and -1 means inverse relationship.

Computing:
Now to compute the convolution/correlation, we use matrix multiplication as a tool, just like  you use variable multiplication for calculating variance or mean in 1 signal case. In multi dimension case(x,y,...) you use matrix usually instead of variables to make your computations easy(especially on computers).

Now in frequency domain:

For calculating, convolution, in time domain, you had to flip and then slide(slide means that each and every element of variable needs to be slided) which is computationally very expensive. Instead, in frequency domain, the same result as that of convolution is achieved just by multiplication and yes when you signal is in freq. domain, you simply multiply the signal(X)-which is matrix with Signal(Y), which is also a matrix.

So, now you will be able to understand that, Yes convolution is same as matrix multiplication(where matrix X and Y matrix of signal) but ONLY IN FREQUENCY DOMAIN.

I hope it helps.

P.S: Simple mulitplication is done in theory, when you are writing equation, but when you have a real signal , then you sample it and it gives a matrix, which you further multiply in matlab.


To answer your 1st question, dilation will enlargen the white spots in your image over the darker spots, which may be useful if you only have tiny dark salting you would like to remove in your image.

However, you are referring to salt noise, so if the spots were also white crumbs in a dark image, then all you are doing is only making the crumbs bigger, which is what you don't want.

What I would recommend is to use either Opening or Closing, (or both) depending on your application. These are a combination of dilations and erosions so that the minor parts of the image are removed while keeping the original portions of the image intact.

For example, here is an image for a penguin (with salting):

Notice how the image has both black and white salting.

Running both an opening and closing algorithm (OpenCV) on this penguin, we remove the black and white salts, and get a nice, smoother image without the spots:

Also, here is my OpenCV code that does this:

#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include <iostream>

using namespace cv;
using namespace std;

int main( )
{
     Mat image,dst;
     image = imread("saltnoise.png", CV_LOAD_IMAGE_GRAYSCALE);

     // Create a structuring element
     int erosion_size;
     cout << "What element size?" << '\n';
     cin >> erosion_size;

     imshow("Salted", image);

     Mat element = getStructuringElement(cv::MORPH_ELLIPSE,
            cv::Size(2 * erosion_size + 1, 2 * erosion_size + 1),
            cv::Point(erosion_size, erosion_size) );

    //Opening and Closing
    erode(image,image,element);
    dilate(image,image,element);

    dilate(image,image,element);
    erode(image,image,element);

     imshow("Unsalted", image);

     waitKey(0);
     return 0;
}


Intuitive understanding of scale-space extrema detection

LoG and DoG (an approximation of LoG) masks can serve as blob detectors. A blob can exist in an image at a number of locations $(x,y)$-coordinates and scales (some parameter; $t$). In some situation where scale space is divided into 3 discrete 'slices' and there are only 'small,' 'medium' and 'large' sized blobs, a 'medium' sized blob will have some response to both the 'small' and 'large' sized detectors. However, a detector whose scale is matched to the scale of the blob will have the largest response.  For this reason, we want to localize the blob by determining the maximum response in both spatial and scale coordinates.


Why DFT is used for approximating CTFT when you can approximate CTFT-integral itself?

Assuming that future vistors won't take the time to read all the comments, I'd like to give a very simple and straightforward interpretation of the discrete Fourier transform (DFT) as an approximation of the continuous-time Fourier transform (CTFT).
First, we need to assume that the signal to be transformed decays sufficiently for large value of $|t|$. (If it doesn't, we can still use the DFT to approximate the CTFT, but the truncation error will be large.) Furthermore, we choose the origin $t=0$ such that the relevant part of the signal appears for $t\ge 0$. Then we can approximate the Fourier integral in the following way:
$$X(j\Omega)=\int_{-\infty}^{\infty}x(t)e^{-j\Omega t}dt\approx \int_{0}^{T}x(t)e^{-j\Omega t}dt\tag{1}$$
with some sufficiently large $T$. The approximation $(1)$ generally results in a truncation error, which can be made small by an appropriate choice of $T$.
Now we can approximate the right-hand side of $(1)$ by a Riemann sum:
$$\int_{0}^{T}x(t)e^{-j\Omega t}dt\approx\sum_{n=0}^{N-1}x(n\Delta t)e^{-j\Omega n\Delta t}\Delta t\tag{2}$$
where $\Delta t$ is the sampling interval and $T=N\Delta t$. The approximation $(2)$ introduces an aliasing error, which can be made small by choosing $\Delta t$ sufficiently small.
Computing the right-hand side of $(2)$ at discrete frequencies $\Omega_k=2\pi k f_s/N$, $k=0,1,\ldots,N-1$, where $f_s=1/\Delta t$ is the sampling frequency,
finally results in
$$\begin{align}X(j\Omega_k)&\approx\sum_{n=0}^{N-1}x(n\Delta t)e^{-j\Omega_k n\Delta t}\Delta t\\&=\frac{1}{f_s}\sum_{n=0}^{N-1}x(n\Delta t)e^{-j\frac{2\pi}{N}nk}\\&=\frac{1}{f_s}\textrm{DFT}\big\{x_d[n]\big\},\qquad x_d[n]=x(n\Delta t)\tag{3}\end{align}$$
Eq. $(3)$ says that the CTFT of $x(t)$ can be approximated on a set of discrete frequencies by the DFT of a sampled version of $x(t)$. Two types of error have to be taken into account: 1. the truncation error, from truncating the signal, and 2. the aliasing error, from sampling the signal.
Note that we can never avoid both errors, because for finite length signals (with zero truncation error) we always get an aliasing error due to the infinite bandwidth of finite length signals, and for perfectly band-limited signals (with zero aliasing error), we always have a truncation error because they can't have finite length. In practice we always have to deal with both types of error.
In sum, answering the question in your question title: under the aforementioned conditions, the DFT is an approximation of the CTFT integral, and it has the advantage that there exist very efficient algorithms (FFT) for its computation.


How can I find expansion coefficients of the y vector in a given basis?

You need to find projection of $\vec{y}$ along unit vectors in the direction of each of the basis vectors $\mathbf v^{(i)}$.

For finding unit vector in the direction of the vector, you just divide the vector by its magnitude.

And, for finding projection along a unit-vector, you just take the dot-product with the unit-vector.

So, the combined steps boil down to doing the following: $$<\vec{y}, \frac{\mathbf v^{i}}{||\mathbf v^{i}||^{2}_{2}}>$$

In your case the basis-vectors are already normalized, meaning they are already unit vectors, so just take the dot-product(inner product) of vector $\vec{y}$ with each of the basis vectors $\mathbf v^i$


Hello again guys after a lot of thought, I finally broke this down and I would like to share it with you hoping it'll help.

So what I did was:

1- I used matlab's pwelch to get the power spectrum like this

[AutoPowerSpectrum ,F] = pwelch(signal, Nsamps, 0, Nsamps, fs, 'power');

you can read about the different parameters for pwelch in matlab's documentation right HERE.

2- I calculated Vpeak squared, because as far as I know pwelch produces Vrms values.

AutoPowerSpectrum = 2.*AutoPowerSpectrum;

3- got the Vpeak by squaring the above equation

AutoPowerSpectrum = sqrt(AutoPowerSpectrum);

4- converted it to dBFs

p_dBFs = 20*(log10(AutoPowerSpectrum./(vPerLevel*2^(N-1))));

where: vPerLevel is the voltage corresponding to each level in my ADC. It is calculated like this Vref/2^N.

and N is the resolution of the ADC.

Finally used plot() to display it. This solution worked with me and gave me the signal in dBFS, i'm opened to edits and suggestions.

Hope this helped.

Cheers.


Denoise Image with Gaussian Noise Using MATLAB / Octave

You're trying to solve what's called Perona Malik Non Linear Diffusion Problem (Sometimes people call it, by mistake, Anisotropic Diffusion).

Anyhow, the simplest code for that is Anisotropic Diffusion (Perona & Malik) on The MATLAB File Exchange.

There is a more advanced (Anisotropic for real) algorithm in Fast Anisotropic Curvature Preserving Smoothing (Also on the File Exchange - Fast Anisotropic Curvature Preserving Smoothing).

Enjoy...


Padding Images To Solve Inverse Filtering Problem?

The thing to understand here is exactly what all your transforms and data fiddling do.

We'll consider a 1D discrete-time signal which is rather easier to explain and think about. The arguments extend easily to higher dimensions.

By taking the DFT of that signal, which has been windowed to a finite number of samples, N, to create a frame, one implicitly assumes that the N-length windowed signal is repeated infinitely in the time-domain, such that the zeroth sample is adjacent to the (N-1)th sample.

The effect of this tiling on what happens when you perform any filtering is two-fold:

1) If there is a big step between the 0th (which is also the Nth) sample and the (N-1)th (which is also the -1th) sample, this is going to introduce lots of high frequency components into the signal which probably don't correspond to anything physical.

2) Filtering by windowing in the Fourier domain is equivalent to convolving with the inverse DFT of that window in the time-domain. Since there is an implicit repeating of the signal in the time domain, the lower samples of the frame "leak" energy into the upper samples of the frame, and vice-versa (alternatively: the ends of the frame are wrapped around and smeared together).

You can't escape these effects, it's the nature of doing discrete Fourier domain filtering. The question is, can one do anything sensible to alleviate them?

The answer depends on knowledge about the signal. You mention zero sampling. Simple zero padding (with no edge roll-off) has the effect of artifically introducing a big edge discontinuity. It's equivalent to convolving the Fourier domain signal with a Sinc function, which is likely going to complicate the required inverse filter. On the other hand, it may alleviate the second of the issues above because now there is a gap between the low samples and the wrapped-around high samples.

If your signal naturally decays to zero at the edge of the frame, then zero padding is sensible and generally safe to separate the tiled versions of the frame.

There are various alternative techniques to simply adding zero samples. One might be to time reverse the frame and concatenate it to the original, making a double length frame. This then guarantees that the signal has no big discontinuities. You can chop the frame in half at the end of the processing to get back what you actually need.

Another technique is to roll off the edges of the time domain signal to zero, typically with something like a raised cosine. This is far less harsh than than simply adding zeros to the time domain. Once the signal has reached zero, you can merrily add more zeros to pad it out if you wish. The length of the roll off is a design parameter, but a good starting point is something like the length of the blurring filter.

It's a case of understanding the signal and the effect of a given tweak as to whether it improves the result or not. I would generally expect to have to do something to alleviate the so called edge effects that you're ostensibly trying to do with zero padding, but I'd try to be a bit more sophisticated about it.

One quick point, of which I'm sure you're aware (so for other people reading this). Simply dividing through by H without any form of regularisation is generally a "bad" thing to do. You need to make sure the signal doesn't become dominated by the very small values of H. That's another post though!


As people who responded in DSPrelated to your question, I believe the fluctuation which are seen are MATLAB plot artifacts not a real envelope. Since you doubt if Fourier analysis could help us to see the spurious signal, I perform a more elaborate spectrum sensing method (MUSIC) to identify if there is other frequency contents in the signal through the code below:

s = 1200*10^6;%Sampling Frequency
T = 1/Fs; % Clock Period
N = 2^12; % No. of samples
t = (0:N-1)*T; % Linear time interval
F = 200.01*10^6; % Required Frequency
phase = 2*pi*F*t; % Consecutive Phase values in rad
sig = sin(phase); % Time Domain Signal
%% SPECTRUM ESTIMATION
X=corrmtx(sig,200,'mod');    % MUSIC SPECTRUM ESTIMATION
[W,P] = rootmusic(X,4);
f=Fs*W./(2*pi); %Frequencies
unique(abs(f)) %Show Frequencies
P              %Power of those frequencies

if you check out the P for each frequency component, you'll see only the main frequency (here 200.01MHz) has non-zero value.


Transfer function of a 1st order HP filter is:

$$H(s)=\frac s{s+wc}=\frac{{\displaystyle\frac s{wc}}}{1+{\displaystyle\frac s{wc}}}\;,where\;wc=2\mathrm\pi\ast\mathrm{fc}$$

&

a = [1 w0];
b = [w0 1e-48]; % ! two coefficients


how to extract meaning from signal (distance estimation)?

In general, you can apply a filter with some time constant T which will be giving you a more reliable RSSI value but even in this case, the RSSI will only be an accurate indication of distance under certain assumptions.

There are two types of filters you could use here, a moving average low pass filter which would effectively produce a smoothed out version of the RSSI waveform or a non-linear modal-mean or median filter. An example of an easily implemented moving average Finite Impulse Response (FIR) filter would be h=ones(N)./N where N is the order of the filter that is directly proportional to the Group Delay (or time constant) of the filter (for instance a [0.25,0.25,0.25,0.25] N=4 filter). Applied over a sliding window of N samples, this filter would return an average, slow changing, estimate of the RSSI curve which you could then associate with distance. Similarly, the non-linear filter is applied over a sliding window of size N but in this case you just estimate the modal mean (i.e. the most frequently occurring RSSI value of the N available in the sample) or the median value. The median filter would tend to remove the spikes and could possibly be a better alternative for this application.

The relationship between RSSI and distance would be established through the Free Space Path Loss but the main assumptions here would be that a) there is a direct line-of-sight between transmitter and receiver and b) the RSSI is accurately estimated by the device. In other words, if you start moving the device in space, it is not necessary that the distance estimation will remain accurate. You might also want to check this link as well.

Hope this helps.


This answer actually provides free code to perform the focus detection.
https://stackoverflow.com/a/3795151/600736


Magnitude response of mirrored (with respect to unit circle) poles and zeros

You're right, the contribution of a zero $z_0=re^{j\phi}$ to the squared magnitude response is
$$\big|e^{j\omega}-re^{j\phi}\big|^2=1-2r\cos(\omega-\phi)+r^2\tag{1}$$
From $(1)$ it is clear that another zero resulting in the same frequency dependence as in $(1)$ must have the same phase angle $\phi$. Assuming $z_1=Re^{j\phi}$ we get
$$\big|e^{j\omega}-Re^{j\phi}\big|^2=1-2R\cos(\omega-\phi)+R^2=R^2\left(1-\frac{2}{R}\cos(\omega-\phi)+\frac{1}{R^2}\right)\tag{2}$$
Comparing $(1)$ and $(2)$ we get $R=1/r$ and a squared magnitude scaling of $R^2=1/r^2$. Consequently, the other zero $z_1$ is given by
$$z_1=\frac{e^{j\phi}}{r}=\frac{1}{z_0^*}\tag{3}$$
and
$$\big|e^{j\omega}-z_0\big|=r\left|e^{j\omega}-\frac{1}{z_0^*}\right|\tag{4}$$


How Do I Measure the Time Duration of a Finite Length Discrete Sequence?

The following graphic should clear up why a useful answer is 5 seconds, consistent with conversions between the discrete-time and continuous-time domains. A good example of this are discrete time approximations to continuous time integration or differentiation where the inclusion of $\Delta T$ (sample period) is made. With all such approximations, the time duration attributed to each sample is one sample interval. For the Forward Euler approximation specifically, each sample covers the time duration up to and not including the next sample. Such as in the first graphic the time duration for a sample at $t=0$ is one sample interval of one second; up to and not including the next sample at $t=1$ when sampled at one sample per second, with a total time duration for the one sample as one second.
The bottom figure then shows how this would apply to five samples starting at $t=0$ resulting in a total time duration of 5 seconds.

The plot below shows the Forward Euler discrete-time approximation for the integration of a continuous time waveform $x(t)$, which would be given as follows in the formula below and demonstrated in the graphic, best showing how each sample carries the time duration of the sample period (and as David has pointed out in the comments, this equally represents the zero-order hold approach in the output of a typical Digital to Analog converter):
$$\int_0^{NT}x(t)dt \approx \sum_{n=0}^{N-1}x[n]T$$

The actual range of time used for each sample is specific to the analog-digital translation approach used, but would cover one sample duration in all cases, leading to a total duration of 5 seconds for the OP's case.


In simple terms, Image shown here speaks for itself.

Before speaking about Fourier Transform black magic lets understand idea behind it. Work of the Mathematician Joseph Fourier demonstrated any arbitrary periodic (this is important) signal can be decomposed into bunch of sine-waves at with their corresponding amplitude and relative phases.

Essentially, he meant basic component to build any signal is a sine-wave. In other words you can build any arbitrary signal exactly if you have information about how many different frequencies are present, amplitude of the each sine wave, and there relative phases. This job can be done by inverse Fourier transform. And as should be clear by now, Fourier Transform does the reverse give it composed time domain signal and it will tell you sine-waves amplitude for corresponding frequencies and there relative phases. Hence, the output is a complex. Which makes sense in polar co-ordinates, which gives you amplitude and phase.

E.G. When you strum a guitar strings, it vibrates at multiple frequencies with different amplitude which makes up a beautiful sound, otherwise you know how a single frequency sounds(watch this :https://www.youtube.com/watch?v=qNf9nzvnd1k). one with the highest amplitude is known as tone.

So, as you can see in the figure attached here, Fourier transform tells you exactly the same information.

  Suppose I have a signal $x:\Bbb R\to \Bbb R$ defined by the formula
  $x(t) = \cos(6\pi x)e^{-\pi x^2}$.

For to obtain Fourier transform of the function you specified it should be function of time (t). So, I will just re-write your function above slightly differently replace x with t, since x(t).

$x(t) = \cos(6\pi t)e^{-\pi t^2}$.

Also, Take a look at this,
http://www.thefouriertransform.com/pairs/rightSidedSinusoids.php

Hope this helps.


How can i define the vector to variable scalar?

Use pointwise vector operations:

% given A,B,C,D,E
A = rand(1,1155);
B = rand(1,1155);
C = rand(1,1155);
D = rand(1,1155);
E = rand(1,1155);
% compute F
F = (A+B)./(C + D.*E);


Calculate Impulse Response of a System From Its Input/Output Waveform via FFT

It would be best to upload the exact data for $x[n]$ and $y[n]$, but the reason I think of, is that either at least one of your signals is not a perfect baseband signal and has considerable amount of energy everywhere, or (this one is more probable) the signals are baseband but your sampling rate is not high enough, hence the Nyquist criterion is not satisfied and aliasing has occurred. In the latter case, try increasing the number of samples per each signal. Note that this is also an approximate, because the DFT package of MATLAB operates on discrete signals, yet with truncation in both time and frequency.


Finding the right measure to compare sound signals in the frequency domain

[EDIT: a related discussion on (power) spectrum centroids can be found in How to calculate the mean/center frequency of the spectrum?]

Let us assume you have a relatively concentrated peak, with noise. I assume your maximum frequency is ok. If too high, you can trim the spectrum to lower frequency, which could help.

The code below gives a simple model, estimates a peak location, and a gross measure of dispersion (or uncertainty). It is based on a weighted average or center of mass. If the peak is narrow, the center of mass (and the median) will be drifted away toward the center of the frequency axis ($f_i$'s). Because even if low, the noise spreads. This is an illustration of Archimedes' lever principle:

  Give me a place to stand on, and I can move the earth.

To achieve your goal, you can give more weight to your spectrum of amplitude $a_i$, and compute an estimated peak location $\bar{f}$:
$$\bar{f} = \frac{\sum_i w_i f_i}{\sum_i w_i}\,. $$

If you choose $w_i = a_i$, you get the standard average. This is not enough. You can increase the peak weight with a power: $w_i = a_i^p$. This is the aim of weightPower = 4 in the code. The location is given by the red o, and the x are the left and right bounds derived from a weighted standard deviation. The higher the $p$, the sharper the estimation, with a risk of major errors if the data differs from the model.

Please report your tests, to adjust the answer. And the result follows. With $p=1$, you get $360$, with $p=2$, you get $243$, with  $p=4$, $\bar{f}=201$, quite close to $200$ in my model. If you need much higher robustness, you can go for weighed medians instead of weighed means.

freqAxis = linspace(0,1000,1000)';
freqPeak = 200;
freqPeakAmplitude = 1;
freqGGDExponent = 2;
freqGGDWidth = 50;
freqNoise = 0.2;

dataSpectrum = exp(-((freqAxis-freqPeak)/freqGGDWidth).^freqGGDExponent)+freqNoise*rand(size(freqAxis));

weightPower = 4;
freqEstLocation = sum((dataSpectrum.^weightPower).*freqAxis)/sum(dataSpectrum.^weightPower);
freqEstDispersion = sqrt(sum((dataSpectrum.^weightPower).*(freqAxis-freqEstLocation).^2)/sum(dataSpectrum.^weightPower));

freqEstAmplitude = interp1(freqAxis,dataSpectrum,freqEstLocation);
clf;hold on
plot(freqAxis,dataSpectrum);
plot(freqEstLocation,freqEstAmplitude,'or');
plot([freqEstLocation-freqEstDispersion freqEstLocation+freqEstDispersion],freqEstAmplitude,'xr');
axis tight;grid on


Is it possible to have high energy compaction in audio signal?

Depends on what you mean exactly by "energy compaction".

Most audio signals are "pink" in nature, so doing any type of frequency transformation will "compact" the energy towards the lower frequencies and most audio codecs use this property. Not sure you get "high" energy compaction: a typical FLAC compression rates are about 30% size reduction or so.

It's more difficult to compact a signal in time through some sort of linear transformation or invertible transformation. This would really depend on what exactly you are trying to do and what the general properties of your signal are.


How can I "smooth" square wave/pulse data to remove/reduce harmonics?

If you are interested in finding the periods, have a look at autocorrelation lag.
Maybe try some xcorr options too.
MATLAB example code:

nh = 47;  % 47 hours.
% Pulse period of 1 hr and make sequence length not a multiple of period
n = 3600*nh + 273; % 1 sec sampling rate
x = zeros(1,n);
% Create a pulse train with varying heights and some noise
x(777:3600:end) = 8 + 2*rand(1,nh);
x = x + rand(1,n);
% Look at lags up to 4 hours
xcorr(x, 4*3600, 'unbiased');

Zero lag is in the middle of the figure.

Note that the consistent heights show reasonably consistent correlation at integer multiples of the base period.
If you want to tolerate timing jitter in the period, you can filter the data with a lowpass FIR filter first, which will smear the peaks in the data and, of course, smear the peaks in the autocorrelation., e.g.,

b = hamming(15);  % Smooth, simple, positive, symmetric, integer delay
x = filter(b,1,x);


FIR definition:

$$ y[n] = \sum_{k=0}^{N} { b_k x[n-k] }  $$

Sinusoid signal definition:

$$ x[n] = M \cos( \alpha n + \phi ) $$

A whole bunch of math:

$$ y[n] = M \sum_{k=0}^{N} { b_k \cos( \alpha (n-k) + \phi ) }  $$

$$ y[n] = M \sum_{k=0}^{N} { b_k [  \cos( \alpha n + \phi ) \cos( \alpha k ) + \sin( \alpha n + \phi ) \sin( \alpha k ) ] }  $$

$$ y[n] = M \cos( \alpha n + \phi ) \sum_{k=0}^{N} { b_k  \cos( \alpha k ) } + M \sin( \alpha n + \phi ) \sum_{k=0}^{N} { b_k \sin( \alpha k ) }  $$

$$ y[n] = A \cos( \alpha n + \phi ) + B  \sin( \alpha n + \phi )  $$

$$ A = M  \sum_{k=0}^{N} { b_k  \cos( \alpha k ) } $$

$$ B = M  \sum_{k=0}^{N} { b_k  \sin( \alpha k ) } $$

$$ y[n] = M_2 \cos( \alpha n + \phi + \theta ) $$

$$ y[n] = M_2 \cos( \theta ) \cos( \alpha n + \phi ) - M_2 \sin( \theta ) \sin( \alpha n  + \phi ) $$

$$ A = M_2 \cos( \theta ) $$

$$ B = -M_2 \sin( \theta ) $$

$$ A^2 + B^2 = M_2^2 = M^2 \left[ \left( \sum_{k=0}^{N} { b_k \cos( \alpha k ) } \right)^2 + \left( \sum_{k=0}^{N} { b_k \sin( \alpha k ) } \right)^2 \right] $$

Your desired equation:

$$ \frac{M}{M_2} = \frac{ 1 }{ \sqrt{ \left( \sum_{k=0}^{N} { b_k \cos( \alpha k ) } \right)^2 + \left( \sum_{k=0}^{N} { b_k \sin( \alpha k ) } \right)^2 } } $$

I think I've done the math right.  The $b_k$s are your FIR coefficients and $\alpha$ is your frequency in radians per sample.  $M_2$ is the amplitude of your smoothed sinusoid and $M$ is the original amplitude.

So you want multiply your smoothed results by $ \frac{M}{M_2} $.

Hope this helps.  I just did this and haven't tested it.

Ced


Spline-Based Adaptive Interpolation Filters?

Look up "Farrow filter". It's essential a set of piecewise polynomial interpolators for the coefficients of an even longer (better) FIR interpolation filter (and without the need for a huge polyphase table).  IIRC, it can be implemented in a straightforward arithmetic hardware pipeline (thus is likely quite suitable for an FPGA).


How would I represent an image using a basis $A$, given a sparse or compressible image?

The textbook algorithm is to minimize

$$\|Ax-m\|_1+\alpha\|x\|_0$$

where $\|x\|_0$ counts the number of non-zero coefficients. Since that is a bad function for the determination of a descent direction, one softens the counting norm to

$$\|x\|_ε=\sum_k|x_k|^ε$$

for some small value of $ε$. For instance, for $|x_k|>\exp(-\tfrac1{\sqrtε})$, which is a very small threshold, one gets $|x_k|^ε>\exp(-\sqrtε)\ge 1-\sqrtε$, so that all not too small and not too large values get mapped to $1$ and only really small values get mapped to zero.

To get a really smooth function, mostly for theoretical purposes, one can smoothen the kink at $0$ by using

$$\|x\|'_ε=\sum_k(\delta^2+|x_k|^2)^{ε/2}$$

and obviously choosing $δ$ small enough so that $x_k=0$ still gets mapped close to $0$, for instance by choosing $δ=\exp(-\tfrac1{\sqrtε^3})$.

Now that this answer is completely orthogonal to the question in its current form, I hope that someone else has more insight into the topic of matrix selection. What I know is this:

In the beginning, before compressed or compressive sensing existed as word, $A$ was a matrix of wavelet or wavelet-like coefficients. The whole idea of wavelet based compression (and JPEG) is to zero out small coefficients (and truncate the binary representation of non-null coefficients). This can be improved upon in one direction by observing and utilizing correlations in the detail coefficients over multiple scales, or in the other direction by making the basis vectors even more uncorrelated and thus reducing those fractal features.

There may also be algorithms that start with a random matrix and modify or reject it based on certain frame conditions.


Why is one of the modulated frequency spikes smaller?

The reason is because of spectral leakage. Your frequencies do not fall exactly on the FFT bins- this causes the energy to appear in nearby adjacent FFT bins as sidelobes. The most spectral leakage occurs when the frequency of a signal falls exactly between two FFT bins. Another way to think of it, is that the signal is not periodic w.r.t. the length of the FFT you are using.

Your two frequencies are appearing at two different positions with respect to the adjacent FFT bins - so they are suffering different amount of leakage. Thus, the peaks have different magnitudes. In the case where the frequencies have the same magnitude, they are located the same distance away from the closest FFT bin.

You can try adjusting your sampling rate and FFT size so that the resulting frequencies will fall exactly on an FFT bin. With the 3000 Hz sampling rate use a 3000 point FFT - corresponding to FFT bins every 1 Hz. With the 4000 Hz sampling rate use a 4000 point FFT - again, corresponding to FFT bins every 1 Hz.

In your 4000 Hz case, you are probably using a 4196 point FFT. The FFT bins are close to being every 1 Hz, so the frequencies at 900 Hz and 1100 Hz are almost exactly on a FFT bin and thus have almost equal amplitude.


Efficient implementation of 2-d circularly symmetric low-pass filter

Approximation by the real part of a weighted sum of separable complex Gaussian component kernels

Figure 1. The proposed scheme illustrated as 1-d real convolutions ($*$) and additions ($+$), for cut-off frequency $\omega_c = \pi/4$ and kernel width $N=41$. Each the upper and the lower half of the diagram is equivalent to taking the real part of a 1-d horizontal and a 1-d vertical pass separable 2-d complex convolution. Except for truncation/cropping, the two final components summed and the resulting kernel are guaranteed circularly symmetrical due to use of complex-number-weighted complex-variance Gaussian 1-d kernels.

I have suggested this approach for this problem, but not tried it before. A discrete-space 2-d convolution kernel (impulse response) $h[x, y]$ that is circularly symmetric (see Eq. 4) about $[0, 0]$, can be approximated to arbitrary accuracy by the real part of a separable (see Eq. 3) complex number weighted sum of Gaussian functions of complex variance:

$$h[x, y] \approx \tilde h[x, y],\tag{1}$$

$$\begin{align}\tilde h[x, y] = &\sum_{m=0}^{M-1}\operatorname{Re}\left(a_m^2e^{b_m(x^2+y^2)}\right)\tag{2}\\
= &\sum_{m=0}^{M-1}\operatorname{Re}\left(a_m e^{b_m x^2}\times a_m e^{b_m y^2}\right)\tag{3}\\
= &\sum_{m=0}^{M-1}\operatorname{Re}\left(a_m^2 e^{b_m r^2}\right)\tag{4}\\
= &\sum_{m=0}^{M-1}e^{\operatorname{Re}(b_m)r^2}\Big({\big(\operatorname{Re}(a_m)^2 - \operatorname{Im}(a_m)^2\big)\cos\big(\operatorname{Im}(b_m)r^2\big)\\
- 2\operatorname{Re}(a_m)\operatorname{Im}(a_m)\sin\big(\operatorname{Im}(b_m)r^2\big)\Big),}\tag{5}\end{align}$$

where $x$ and $y$ are integer horizontal and vertical coordinates, $r = \sqrt{x^2 + y^2}$ is the radius or distance between $[0, 0]$ and $[x, y]$, $M$ is the approximation order, and $a_m$ and $b_m$ are complex amplitude and variance related coefficients.

Optimized approximation, $\omega_c = \pi$

Eq. 5 with $r^2 = x^2 + y^2$ can be numerically optimized by global optimization to approximate a given $h[x, y]$ in spatial domain within the square $-(N-1)/2 \le x \le (N-1)/2$, $-(N-1)/2 \le y \le (N-1)/2$ with odd width $N$ of the square kernel. It suffices to do the optimization over a triangle $0 \le y \le (N-1)/2$, $y \le x \le (N-1)/2$. If the cost function is a sum of squares of errors or another sum of errors at each $(x, y)$, then cumulative $2\times$ error multipliers should be applied for each satisfied condition $x > 0$, $y > 0$, $x > y$.

The following C++ program optimizes the approximate kernel $\tilde h$ in spatial domain least squares sense using Differential Evolution from an optimization library. The target kernel $h$ is windowed using a 2-d rotated cosine window. The source code includes pre-optimized approximations for cut-off frequency $\omega_c = \pi$, kernel size $N=41$, and approximation orders $1 \le m \le 5$.

// -*- compile-command: "g++ -Wno-unused-result -march=native -O3 -ffast-math optitest.cpp opti.cpp" -*-
// This file is optitest.cpp
#include <stdio.h>
#define _USE_MATH_DEFINES
#include <math.h>
#include "opti.hpp"
#include "keyboard.h"

double optimal0[4] = { // Dummy starting point for optimization
  0,0,0,0
};

double optimal1[1*4] = {//bestcost=0.03593696209385285195 Least squares
  0.90954095944946312,0.21237789645057972,-0.72663606661376057,0.76147736559107793
};

double optimal2[2*4] = {//bestcost=0.01072419992771576382 Least squares
  -0.14935266179235707,0.28729650415029556,-0.08227954155725942,0.36475170816661134,1.03225137134955114,0.47158356759095016,-0.60039465413238513,0.47344404338750434
};

double optimal3[3*4] = {//bestcost=0.00322191268986909149 Least squares
  -0.00646998371362690,-0.15983965263134517,-0.02321635125852370,0.20289085505437962,-0.30298121364781033,0.34395623806288950,-0.10323939264492392,0.32564190139614663,1.13148495464950427,0.66639187966234981,-0.54567361845465989,0.37519342198291905
};

double optimal4[4*4] = {//bestcost=0.00089888657916611925 Least squares
  0.04481214216545543,0.08998210660870602,-0.01130470880820374,0.14003687599487286,0.02506113778123746,-0.24720956687380186,-0.03419309947921907,0.19643391040650712,1.27604032320787875,0.91159589755622883,-0.49585695598458995,0.29480666144138823,0.48259638667856242,-0.42280882153371496,-0.12730997070336811,0.29988730534029784
};

double optimal5[5*4] = {//bestcost=0.00017259634849055045 Least squares
  -0.06413980110992069,-0.03272558679644168,-0.00736621171073370,0.10943881870260203,0.08213936888117918,0.16226536153011967,-0.01975126456281387,0.13958603391531316,-0.07270617897425770,0.37473199916953354,-0.04702984945995840,0.18925309227383197,1.43523098843984531,1.17056740170289952,-0.45755771004055446,0.24332707322808175,0.69700596730108921,-0.48854609666142051,-0.14144974773647198,0.27672159791886242
};

const int numOptimalKernels = 5;
double *optimalKernels[numOptimalKernels + 1] = {optimal0, optimal1, optimal2, optimal3, optimal4, optimal5};

class CircularLPProblem : public Opti::Problem {
private:
  double *minimum;
  double *maximum;
  int numComponents;
  int numSamples;
  int *r2s;
  int *multipliers;
  double *target;
public:

  int getNumDimensions() {
    return numComponents*4;
  }

  double *getMin() {
    return minimum;
  }

  double *getMax() {
    return maximum;
  }

  double costFunction(double *params, double compare) {
    for (int m = 0; m < numComponents; m++) {
      params[4*m + 2] = -fabs(params[4*m + 2]);
      params[4*m + 3] = fmod(fabs(params[4*m + 3]), M_PI);
    }
    for (int m = 0; m < numComponents - 1; m++) {
      if (params[4*m + 3] > params[4*(m + 1) + 3]) {
        for (int k = 0; k < 4; k++) {
          double temp = params[4*m + k];
          params[4*m + k] = params[4*(m + 1) + k];
          params[4*(m + 1) + k] = temp;
        }
      }
    }
    double cost = 0;
    for (int k = 0; k < numSamples; k++) {
      double r2 = r2s[k];
      double trial = 0;
      for (int m = 0; m < numComponents; m++) {
        //  trial += exp(params[4*m + 2]*r2)*(params[4*m + 0]*cos(params[4*m + 3]*r2) + params[4*m + 1]*sin(params[4*m + 3]*r2)); // Max absolute error
        trial += exp(params[4*m + 2]*r2)*((params[4*m + 0]*params[4*m + 0] - params[4*m + 1]*params[4*m + 1])*cos(params[4*m + 3]*r2) - 2*params[4*m + 0]*params[4*m + 1]*sin(params[4*m + 3]*r2)); // Least squares
      }
      /*      if (fabs(trial - target[k]) > cost) { // Max absolute error
              cost = fabs(trial - target[k]);
              }*/
      cost += (trial - target[k])*(trial - target[k])*multipliers[k]; // Least squares
      if (cost > compare) {
        return cost;
      }
    }
    return cost;
  }

  // numComponents = number of components
  // N = kernel width
  // omega_c = cutoff frequency (radians)
  // Parameter vector: Re(a_0), Im(a_0), Re(b_0), Im(b_0), Re(a_1), ...
  CircularLPProblem(int numComponents, int N, double omega_c = M_PI): numComponents(numComponents) {
    numSamples = 0;
    for (int y = 0; y < (N-1)/2 + 1; y++) {
      numSamples += (N-1)/2 + 1 - y;
    }
    r2s = new int[numSamples];
    multipliers = new int[numSamples];
    target = new double[numSamples];
    int k = 0;
    for (int y = 0; y < (N-1)/2 + 1; y++) {
      for (int x = y; x < (N-1)/2 + 1; x++) {
        r2s[k] = x*x + y*y;
        target[k] = omega_c*j1(omega_c*sqrt(x*x + y*y))/(2*M_PI*sqrt(x*x + y*y));
        double window = cos(M_PI/2*sqrt(pow(x/((N - 1)/2 + 1.0), 2) + pow(y/((N - 1)/2 + 1.0), 2)));
        if (window < 0) {
          target[k] = 0;
        } else {
          target[k] *= window;
        }
        multipliers[k] = ((x > 0) ? 2 : 1) * ((y > 0) ? 2 : 1) * ((x > y) ? 2 : 1);
        k++;
      }
    }
    target[0] = omega_c*omega_c/(4*M_PI);
    minimum = new double[4*numComponents];
    maximum = new double[4*numComponents];
    k = 0;
    for (int i = 0; i < 4*numComponents; i++) {
      minimum[i] = optimalKernels[numComponents - 1][k]-pow(0.1, numComponents - 1);
      maximum[i] = optimalKernels[numComponents - 1][k]+pow(0.1, numComponents - 1);
      k++;
      if (k >= (numComponents-1)*4) {
        k -= 4;
        if (k < 0) {
          k = 0;
        }
      }
    }
  }

  ~CircularLPProblem() {
    delete[] minimum;
    delete[] maximum;
    delete[] r2s;
    delete[] multipliers;
    delete[] target;
  }
};

int main()
{
  INITKEYBOARD;

  CircularLPProblem problem(1, 41, M_PI); // Parameterize this!

  Opti::Strategy *optimizer;

  //optimizer=new Opti::G3(&problem, 2000);
  optimizer=new Opti::DE(&problem, 2000);

  printf("\nOptimizing...\n\n");

  for(int t = 0;; t++) {
    double bestcost = optimizer->evolve();
    if (!(t % 1000)) {
      printf("gen=%d, bestcost=%.20f, average=%.20f\n", t, bestcost, optimizer->averageCost());
      if (kbhit()) {
        printf("Parameter vector printout:\n");
        problem.print(optimizer->best());
        if (getch() == 27) break;
        getch();
      }
    }
  }
  delete optimizer;

  DEINITKEYBOARD;
  return 0;
}

Figure 2. Sum-of-squares cost of the best found approximate kernel $\tilde h$ as function of approximation order $M$, for $\omega_c = \pi$, $\omega_c = \pi/2$, and $\omega_c = \pi/4$, and $N=41$. A good indicator that these are the globally optimal parameterizations is the steady decrease of error as $M$ is incremented, except for $\omega_c = \pi/4$, $M=3$ for which we might (or not) have a sub-optimal solution.

This Python script generates the "exact" (windowed) kernel for comparison and implements the approximations:

import matplotlib.pyplot as plt
from scipy import special
import numpy as np
import scipy.ndimage
import skimage
import time

def circularLowpassKernel(omega_c, N):  # omega = cutoff frequency in radians (pi is max), N = horizontal size of the kernel, also its vertical size, must be odd.
  kernel = np.fromfunction(lambda x, y: omega_c*special.j1(omega_c*np.sqrt((x - (N - 1)/2)**2 + (y - (N - 1)/2)**2))/(2*np.pi*np.sqrt((x - (N - 1)/2)**2 + (y - (N - 1)/2)**2)), [N, N])
  kernel[(N - 1)//2, (N - 1)//2] = omega_c**2/(4*np.pi)
  return kernel

def rotatedCosineWindow(N):  # N = horizontal size of the targeted kernel, also its vertical size, must be odd.
  return np.fromfunction(lambda y, x: np.maximum(np.cos(np.pi/2*np.sqrt(((x - (N - 1)/2)/((N - 1)/2 + 1))**2 + ((y - (N - 1)/2)/((N - 1)/2 + 1))**2)), 0), [N, N])

N = 41  # Kernel width
M = 1  # Number of components
omega_c = np.pi  # Cutoff frequency <= np.pi
kernel = circularLowpassKernel(omega_c, N)*rotatedCosineWindow(N)

def saveKernel(name, kernel):
  plt.imsave(name+'.png', plt.cm.bwr(plt.Normalize(vmin=-kernel.max(), vmax=kernel.max())(skimage.transform.rescale(kernel, 4, 0))))
  absF = np.abs(np.fft.fftshift(np.fft.fft2(kernel)))
  plt.imsave(name+'_dft.png', plt.cm.Greys(plt.Normalize(vmin=0, vmax=absF.max())(skimage.transform.rescale(absF, 4, 0))))

saveKernel('exactpi', kernel)

plt.imsave('bwrkey.png', plt.cm.bwr(np.repeat([(np.arange(41*4)/(41*4-1))], 16, 0)))
plt.imsave('Greyskey.png', plt.cm.Greys(np.repeat([(np.arange(41*4)/(41*4-1))], 16, 0)))

def reComponentKernel(N, param):
  x = np.arange(N)-(N - 1)//2
  return np.exp(param[2]*x**2)*(param[0]*np.cos(param[3]*x**2) - param[1]*np.sin(param[3]*x**2))

def imComponentKernel(N, param):
  x = np.arange(N)-(N - 1)//2
  return np.exp(param[2]*x**2)*(param[1]*np.cos(param[3]*x**2) + param[0]*np.sin(param[3]*x**2))

optimal1pi = [0.90954095944946312,0.21237789645057972,-0.72663606661376057,0.76147736559107793]
optimal2pi = [-0.14935266179235707,0.28729650415029556,-0.08227954155725942,0.36475170816661134,1.03225137134955114,0.47158356759095016,-0.60039465413238513,0.47344404338750434]
optimal3pi = [-0.00646998371362690,-0.15983965263134517,-0.02321635125852370,0.20289085505437962,-0.30298121364781033,0.34395623806288950,-0.10323939264492392,0.32564190139614663,1.13148495464950427,0.66639187966234981,-0.54567361845465989,0.37519342198291905]
optimal4pi = [0.04481214216545543,0.08998210660870602,-0.01130470880820374,0.14003687599487286,0.02506113778123746,-0.24720956687380186,-0.03419309947921907,0.19643391040650712,1.27604032320787875,0.91159589755622883,-0.49585695598458995,0.29480666144138823,0.48259638667856242,-0.42280882153371496,-0.12730997070336811,0.29988730534029784]
optimal5pi = [-0.06413980110992069,-0.03272558679644168,-0.00736621171073370,0.10943881870260203,0.08213936888117918,0.16226536153011967,-0.01975126456281387,0.13958603391531316,-0.07270617897425770,0.37473199916953354,-0.04702984945995840,0.18925309227383197,1.43523098843984531,1.17056740170289952,-0.45755771004055446,0.24332707322808175,0.69700596730108921,-0.48854609666142051,-0.14144974773647198,0.27672159791886242]

class SeparableCircularLowpassFilter:
  def __init__(self, N, coefs):
    self.N = N
    self.reKernels = []
    self.imKernels = []
    for i in range(len(coefs)//4):
      self.reKernels.append(np.array([reComponentKernel(N, coefs[i*4:])]))
      self.imKernels.append(np.array([imComponentKernel(N, coefs[i*4:])]))
  def filter(self, x):
    reZ = scipy.ndimage.convolve(scipy.ndimage.convolve(x, self.reKernels[0].transpose()), self.reKernels[0]) - scipy.ndimage.convolve(scipy.ndimage.convolve(x, self.imKernels[0].transpose()), self.imKernels[0])
    for i in range(1, len(self.reKernels)):
      reZ += scipy.ndimage.convolve(scipy.ndimage.convolve(x, self.reKernels[i].transpose()), self.reKernels[i]) - scipy.ndimage.convolve(scipy.ndimage.convolve(x, self.imKernels[i].transpose()), self.imKernels[i])
    return reZ

filter1pi = SeparableCircularLowpassFilter(N, optimal1pi) # 1 component kernel
filter2pi = SeparableCircularLowpassFilter(N, optimal2pi) # 2 component kernel
filter3pi = SeparableCircularLowpassFilter(N, optimal3pi) # 3 component kernel
filter4pi = SeparableCircularLowpassFilter(N, optimal4pi) # 4 component kernel
filter5pi = SeparableCircularLowpassFilter(N, optimal5pi) # 5 component kernel

x = np.zeros([N, N])  # Input image
x[N//2, N//2] = 1     # (unit impulse)
# x = plt.imread('sample.tif').astype(float)[:,:,1]  # (green channel of some image from file)

t0 = time.time()
Z = scipy.ndimage.convolve(x, kernel)  # Exact
t1 = time.time()
print(t1-t0)

t0 = time.time()
reZ1pi = filter1pi.filter(x)  # 1 component kernel
t1 = time.time()
print(t1-t0)
saveKernel('reZ1pi', reZ1pi)

t0 = time.time()
reZ2pi = filter2pi.filter(x)  # 2 component kernel
t1 = time.time()
print(t1-t0)
saveKernel('reZ2pi', reZ2pi)

t0 = time.time()
reZ3pi = filter3pi.filter(x)  # 3 component kernel
t1 = time.time()
print(t1-t0)
saveKernel('reZ3pi', reZ3pi)

t0 = time.time()
reZ4pi = filter4pi.filter(x)  # 4 component kernel
t1 = time.time()
print(t1-t0)
saveKernel('reZ4pi', reZ4pi)

t0 = time.time()
reZ5pi = filter5pi.filter(x)  # 5 component kernel
t1 = time.time()
print(t1-t0)
saveKernel('reZ5pi', reZ5pi)

$M=1$
$M=2$
$M=3$
$M=4$
$M=5$

exactFigure 3. Kernels and the absolute value of their discrete Fourier transform (DFT), enlarged by a factor of 4 to make individual pixels visible. Top to bottom: Optimal approximating kernels $\tilde h$ with $M = 1$, $2$, $3$, $4$ and $5$ complex separable components, ideal kernel $h$ for $\omega_c = \pi$ and $N = 41$. Color keys: kernel: blue negative, white zero, red positive (normalized); abs DFT: white zero, black maximum.

Optimized approximation, $\omega_c = \pi/2$

Approximation of lower-cutoff kernels works better with this approach. I think this is because the number of ripples that fit in the windowed target kernel are reduced. For a high-cutoff filter then perhaps the kernel width $N$ could be reduced, as this would give the same number of ripples and frequency-domain approximation accuracy. But that would then favor a direct 2-d implementation (see Fig. 6). The C++ source code in the above was modified (not shown) for $\omega_c = \pi/2$ (and later for $\omega_c = \pi/4$) and the obtained parameters were used in a Python implementations of the approximate isotropic low-pass filters: (continued from the previous script)

omega_c = np.pi/2
kernelpi2 = circularLowpassKernel(omega_c, N)*rotatedCosineWindow(N)
saveKernel('exactpi2', kernelpi2)

optimal1pi2 = [0.44103810622146067,0.08998875769710178,-0.17010258583392401,0.19960767673288432]
optimal2pi2 = [-0.07233719880423649,0.14289630144713414,-0.02205699413927855,0.09165233018125875,0.51013245392078410,0.22632987351129516,-0.15215724343836151,0.12131467270512424]
optimal3pi2 = [-0.03867959516035375,-0.07754823299868645,-0.00905504286350691,0.05573851697983074,0.60460672854618647,0.36747347219951876,-0.13550005930232881,0.08061493799161984,-0.15528642640407436,0.23307152837452039,-0.03561076440525033,0.08259420496099962]
optimal4pi2 = [0.06323093028956613,0.02993129946018375,-0.00937470404526276,0.03826932634049100,0.78321346985185014,0.58823521191007977,-0.11611698422922974,0.05181454568143690,-0.07541120991444331,-0.16653642294319324,-0.01587217217852400,0.05400812761661748,-0.30454011473801174,0.38069299325271988,-0.04600795167685135,0.07149320869042795]

filter1pi2 = SeparableCircularLowpassFilter(N, optimal1pi2) # 1 component kernel
filter2pi2 = SeparableCircularLowpassFilter(N, optimal2pi2) # 2 component kernel
filter3pi2 = SeparableCircularLowpassFilter(N, optimal3pi2) # 3 component kernel
filter4pi2 = SeparableCircularLowpassFilter(N, optimal4pi2) # 4 component kernel

reZ1pi2 = filter1pi2.filter(x)
saveKernel('reZ1pi2', reZ1pi2)
reZ2pi2 = filter2pi2.filter(x)
saveKernel('reZ2pi2', reZ2pi2)
reZ3pi2 = filter3pi2.filter(x)
saveKernel('reZ3pi2', reZ3pi2)
reZ4pi2 = filter4pi2.filter(x)
saveKernel('reZ4pi2', reZ4pi2)

$M=1$
$M=2$
$M=3$
$M=4$

exactFigure 4. Kernels and the absolute value of their discrete Fourier transform (DFT), enlarged by a factor of 4 to make individual pixels visible. Top to bottom: Optimal approximating kernels $\tilde h$ with $M = 1$, $2$, $3$, and $4$ complex separable components, ideal kernel $h$ for $\omega_c = \pi/2$ and $N = 41$. Color keys: kernel: blue negative, white zero, red positive (normalized); abs DFT: white zero, black maximum.

Optimized approximation, $\omega_c = \pi/4$

In Python: (continued from the first Python script)

omega_c = np.pi/4
kernelpi4 = circularLowpassKernel(omega_c, N)*rotatedCosineWindow(N)
saveKernel('exactpi4', kernelpi4)

optimal1pi4 = [0.22210250507512549,0.04738076096442791,-0.04565916667208288,0.04924268271256500]
optimal2pi4 = [0.28445762468654334,0.14050386040493884,-0.03988634664177378,0.02123143861749057,-0.01717623293554685,0.11158938041615812,-0.01191733840600387,0.02614013272136841]
optimal3pi4 = [-0.00332428227606697,0.01722626235358831,-0.00516799783527921,0.00727825629264402,0.30296628929626013,0.16328004086794692,-0.03813290200362767,0.01805103148816284,-0.02296006483369628,0.12768470633404483,-0.01299438181561595,0.02558273151301529]

filter1pi4 = SeparableCircularLowpassFilter(N, optimal1pi4) # 1 component kernel
filter2pi4 = SeparableCircularLowpassFilter(N, optimal2pi4) # 2 component kernel
filter3pi4 = SeparableCircularLowpassFilter(N, optimal3pi4) # 3 component kernel

reZ1pi4 = filter1pi4.filter(x)
saveKernel('reZ1pi4', reZ1pi4)
reZ2pi4 = filter2pi4.filter(x)
saveKernel('reZ2pi4', reZ2pi4)
reZ3pi4 = filter3pi4.filter(x)
saveKernel('reZ3pi4', reZ3pi4)

$M=1$
$M=2$
$M=3$

exactFigure 5. Kernels and the absolute value of their discrete Fourier transform (DFT), enlarged by a factor of 4 to make individual pixels visible. Top to bottom: Optimal approximating kernels $\tilde h$ with $M = 1$, $2$ and $3$ complex separable components, ideal kernel $h$ for $\omega_c = \pi/4$ and $N = 41$. Color keys: kernel: blue negative, white zero, red positive (normalized); abs DFT: white zero, black maximum.

Computational complexity

Complexity as the number of 1-d dot products of length $N$ per pixel is $4M$ for the approximations and $N$ for the naive 2-d convolution implementation. The approximations suffer from an additional performance overhead due to the use of additional memory for storing intermediate results.

Figure 6. Complexity as the number of real multiplications per pixel vs. kernel width $N$. In reality the approximations are not quite as efficient, see Fig. 7. Possible advantages from symmetry of the kernels has not been taken into account.

Figure 7. Actual measured run time in seconds for 2-d convolution of a 5424×3636-pixel single-channel image by the ideal kernel $h$ and by separable approximating kernels $\tilde h$ with $M = 1$, $2$, $3$, $4$ and $5$ complex separable components, for $N = 41$, using Python's scipy.ndimage.convolve.

Further ideas

In order to adjust $\omega_c$ without kernel parameter re-optimization, a nearly continuous approximating kernel could be optimized by using a small $\omega_c$ and a large $N$. The result could then be sampled using interpolation.

It would be interesting to try direct optimization of the complex kernels without a parameterization that enforces circular symmetry. Perhaps the approximation could even be optimized as a sum of $M$ real separable kernels, similar to this Mathematics Stack Exchange question. Also minimization of (least squares) error in the frequency domain is something to try.


Can a digital image be considered as a tensor ? or have a tensor like structure?

An isotropic tensor is one invariant under rotations. An isotropic kernel is also invariant under rotations, but the way the rotation acts on the object is very different.

Tensors transform as linear forms, whereas images are just rotated like a scalar function of coordinates and the coordinates are rotated. So if you come from continuum mechanics, the image/kernel transforms like the pressure field.


Demodulation of BFSK signal with frequency hopping using Goertzel Algorithm

Guys I am really sorry about this stupid question. Just shifting the center frequency of the spectrum from 0Hz to some higher value ( I used 3000Hz as the center frequency) and then calculating the power using goertzel algorithm at the shifted frequencies, solved the problem. I forgot about the 2-sided nature of the frequency spectrum.


Advantage of using FIR Polyphase Decimator over just computing FIR output every D samples

One advantage of the polyphase filter is when hardware resources are limited. The polyphase approach replaces a $L$ tap filter with $N$ filter sets of $L/N$ taps. Once you've done that you can just use a single $L/N$ filter set and swap coefficients.


I'm going to go out on a limb and answer this based on my gut instinct.

Let's suppose $H(z)$ is casual and has a single zero (a) and a single pole (b):

$$
H(z) = \frac{(1-az^{-1})}{(1-bz^{-1})}
$$

if the system is casual, then we know that the sequence is right-sided with ROC that is exterior of the circle. Thus, the ROC of H(z) is:

$$b < |z|$$

Now for the part I was uncertain about: the ROC of the inverse system G(z).

$$
G(z) = \frac{1}{H(z)} = \frac{(1-bz^{-1})}{(1-az^{-1})}
$$

When we take the inverse of H(z), it is necessary for a portion of the ROC of G(z) to overlap with a portion of the ROC of H(z) in order for an inverse function G(z) to exist. The overlapping doesn't necessary need to include the entire ROC just a portion of it.  This is the constraint that determines if G(z) has ROC that is casual / exterior-of-circle,   or anti-casual / interior-of-circle.  Its a matter of noting the poles of G(z) and checking the following cases to determine direction of the inequality for ROC:

(1) does the exterior of the largest pole of G(z) overlap with a portion of ROC of H(z)? if so then this represents a valid causal ROC (that is an exterior circle). ROC: a < |z|

(2) does the interior of the smallest pole of G(z) overlap with a portion of the ROC of H(z)? if so then this represents a valid anti-casual ROC (that is an interior circle). ROC: |z| < a

it is entirely possible that both a causal inverse and an anti-casual inverse both exist and are valid.

Another consideration for choosing a valid inverse function of H(z) is to make sure that the inverse function is stable, ie. the ROC includes the unit circle, which enables the system to be realizable in real hardware, and is more than a mere mathematical abstraction.


I believe this is due to your choice of using a rectangular window, which has the best frequency resolution (1 bin) but worst spectral leakage (dynamic range). The spectral leakage is of the form of a $\textrm{sinc}$ function (with nulls at $1/T$ where $T$ is your window length in time), which has an envelope that rolls off very slowly ($1/f$) in frequency. What I believe you see above is the cumulative effect of the $\textrm{sinc}$ function in frequency due to your window convolving with each tone.

My favorite window choice is the Kaiser window kaiser(N,beta), as you can control the precision and dynamic range trade directly through the use of beta. Beta factors I typically use range from 4 to 12 but you can experiment with it and see what works best for you.

Try this to see the window performance (MATLAB / Octave) in the frequency domain (showing the main lobe width (precision) and sidelobe level (dynamic range).

wind = kaiser(N, beta)
freqz(wind)


Implementing FIR Hilbert Transformer in Python

You can't implement an ideal Hilbert transformer: its impulse response is non causal and infinite in time.
So you can only implement an "approximate" one and the best way to do this depends on the specific requirement of your application: what's your frequency range of interest, how much magnitude and/or phase deviation can you tolerate, are you constrained in latency or computational resources, do you need to maintain causality and transients, etc.
Python has a Hilbert transform function https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.hilbert.html but without understanding your requirements first, it will be difficult to control the parameters properly.
Very good article on Hilbert transform and analytical signals in the discrete domain: http://andrewduncan.net/air/


Here is a visual answer why it is not working. The overlap is not even symmetrical in the last step so the scheme can't work.

Even if you make it symmetrical by making the filter zero-delay, the amount of stuff that gets summed together varies by time, so it can't be right.


Where can I find a table of $z$-domain coefficients for Butterworth filters?

You can't find a table of $z$-domain coefficients that is analogous to the tables of Laplace-domain Butterworth transfer functions because it's impossible. For a second-order Butterworth filter, for example, there is no single corresponding $z$-domain transfer function, in contrast to the continuous-time case.

The reason: discrete-time transfer functions don't cover an infinitely-long frequency axis like a continuous-time system does. When creating a digital approximation to an analog prototype, techniques like the bilinear transform are used to map the frequency response from the $s$-domain transfer function evaluated along the (infinitely long) $j\omega$ axis to the $z$-domain's (finite-length) unit circle. This obviously can't be done exactly, so the various methods for mapping an analog filter to a digital one are approximations.

One factor that plays into the quality of the approximation is the sample rate of the target digital filter. In the paper you linked, you noted that the sampling interval $T$ cancelled out in the calculations of the digital filter coefficients, which is true. However, the sample interval is still required in order to perform the (necessary) frequency prewarping. Recall the prewarping formula:

$$
\omega_p = \frac{2}{T}\tan{\left(\frac{\omega T}{2}\right)}
$$

While the rest of the $T$'s may cancel out after performing the substitution $s = \frac{2}{T}\frac{z-1}{z+1}$, the $T$ in the argument to the tangent function cannot go away. Therefore, the resulting coefficients will always be a function of the sample rate $T$, so your hope of creating a simple table of digital Butterworth filter coefficients invariant of sample rate is stymied.

A couple other notes regarding this:

The prewarping step is essential to ensure that particular frequencies of interest (i.e. the cutoff frequency for a Butterworth filter) are preserved through the approximation process. Without prewarping, then the stretching of the frequency response that occurs when mapping from an infinitely-long frequency axis to the unit circle would cause the cutoff to appear at a different frequency than desired in the digital filter.
The cookbook linked in niaren's answer also demonstrates the dependence upon the sample rate: you have to trudge through the equations a bit, but you'll see factors in the listed coefficients like $\cos(\omega_0)$ and so on (this is really just a way of slightly hiding the prewarping process). The cookbook defines $\omega_0$ as:
$$
\omega_0 = \frac{2\pi f_0}{F_s}
$$
so when it all comes down to it, you'll have to define your sample rate and prewarp any frequencies of interest in order to design a Butterworth filter approximation. It's not particularly difficult, as the cookbook shows; I'm not sure what the motivation for avoiding those steps is.


Is Differentiation as a system, is an invertible system?

Yes, your understanding is correct. The differentiation is not an invertible system. And, the reason is precisely, because of the possible unknown constant which gets nulled by differentiation.

Basically, the differentiation of $f(t)$ and $(f(t) + c)$ is indistinguishable $f^{'}(t)$, which makes the system, many-to-one, and hence non-invertible.


Why discrete wavelet transforms use sampling rate 2 and need signal length to be a power of 2?

I can't comment because of reputation but if i'm not mistaken, the need for DWT to be supplied an input that is length power of 2 is independent of it's subsequent levels being downsampled by a power of two. The input does not have to be a power of two but you'll have to deal with problems at the edges of the signal then.
Here is a good link for reference that addresses padding of a signal for wavelets:
https://paos.colorado.edu/research/wavelets/faq.html#padding
Basically states that to prevent problems at the edges of the tranform, the inputs are usually padded.
I've also seen places mention other decimation schemes but it begs the question as to why they would be necessary?
For every level of decomposition in a DWT, the signal is downsampled by 2 because the subsequent level of decomposition is limited to half the bandwidth of the previous level. Allowing for half of the samples to be discarded.


Hilbert transform to compute signal envelope?

The Hilbert transform is used to calculate the "analytic" signal. See for example http://en.wikipedia.org/wiki/Analytic_signal. If your signal is a sine wave or an modulated sine wave, the magnitude of the analytic signal will indeed look like the envelope. However, the computation of the Hilbert transform is not trivial. Technically it requires a non-causal FIR filter of considerable length so it will require a fair amount of MIPS, memory and latency.

For a broad band signal, it really depends on how you define "envelope" for your specific application. For your application of dynamic range compression you want a metric that is well correlated with the the perception of loudness over time. The Hilbert Transform is not the right tool for that.

A better option would be to apply an A-weighted filter (http://en.wikipedia.org/wiki/A-weighting) and then do a lossy peak or lossy RMS detector. This will correlate fairly well with perceived loudness over time and is relatively cheap to do.


Alright, after some days of staring at this problem I hope I can provide a bit of guidance to the next poor soul.

Yes the scaling is different. Compared to scipy.fftpack.dct the DC term is $\frac{1}{2}$ and the other terms $\frac{\sqrt{2}}{2}$. But apparently it all nicely cancels out in the inverse transform.
The inverse input order is exactly as they come out: bit reversed. Literally as if you flip around the figure and connect the lines. And yes, you are correct that the sin gets a minus. I'm seeing a scaling factor of 8, FWIW.
The output order of the inverse is the same as the input order of the forward transform. So $x[n]=\frac{\mathrm{IDCT}(\mathrm{DCT}(x[n]))}{8}$

Also note that there is an error in the graph and it's $\sqrt{2}c6$ for the even-side rotation block.


Signal reconstructed using IFFT: amplitude scale larger than input

It's because you are probably adding back overlapping portions together, which adds same signal $N$ times where it seems $N \approx 10$ in your case. To test your algorithm use a zero overlapping case and see that IFFT produces the signal without a gain factor.


Why is the state space representation, $dx/dt = A x(t) + B u(t)$ linear?

What you have listed is a state equation.  Usually a state-space representation has two equations:

(i) $\dot{x}(t)=Ax(t)+Bu(t)$
(ii) $y(t)=Cx(t)+Du(t)$

One takes the input $u(t)$ and generates the state $x(t)$ and another that takes the input $u(t)$ and state $x(t)$ and generates the output $y(t)$.  You are missing the second equation.

That aside, the equation that you have is basically a 1st order linear differential equation.  It's important to note, here, that the input is not $x(t)$, it is $u(t)$.  $x(t)$ is the state variable.

I will examine just the first half (input to state) and leave the second half (input and state to output) to you.  I will remark, though, that if $x(t)$ is linear with respect to $u(t)$ then $y$ in equation (ii) above necessarily has to be linear by virtual of adding two linear terms together, so the second half is really trivial.

The first half is:

(1) $\dot{x}(t)=Ax(t)+Bu(t)$

Can you show that it is linear?

I will show additivity and leave homogeneity to you.  Say you have inputs $u_1$, $u_2$, and $u_3=u_1+u_2$.

One way to do this is to solve (1).  As it is a 1st order ODE, it is not too difficult, but let's be lazy.

Let $x_1$ be the state generated by input $u_1$ and $x_2$ be the state generated from input $u_2$, then we have

(2) $\dot{x_1}(t)=Ax_1(t)+Bu_1(t)$

and

(3) $\dot{x_2}(t)=Ax_2(t)+Bu_2(t)$

adding (2) and (3) gives

(4) $\dot{x_1}(t)+\dot{x_2}(t)=Ax_1(t)+Ax_2(t)+Bu_1(t)+Bu_2(t)$

From linearity of derivative, we get that:

(5) $\dot{x_3}(t)=A(x_3(t))+B(u_3(t))$

where $x_3 = x_1+x_2$.

So input $u_3 = u_1+u_2$ generates state variable $x_3 = x_1+x_2$, showing that additivity holds for the state equation.


The error between the theoretical value and the simulation value of imaginary part of the FFT results

Judging from your question, I think that you're confusing two related but different Hilbert relations. The first one holds for causal signals, for which the real and imaginary parts of their Fourier transforms are related via the Hilbert transform. The other one holds for analytic signals (whose spectrum is zero for negative frequencies), for which the real and imaginary parts of the time domain signal are related via the Hilbert transform.

So if you have a real valued function, as described in your question, and you take its Fourier transform, then the imaginary part of the Fourier transform can only be obtained from its real part via the Hilbert transform, if the signal is causal. Also note that the corresponding time domain signal is generally not analytic (because in that case its spectrum must be zero for negative frequencies and, consequently, the time domain signal must be complex-valued).


python; audio; detecting silence in audio signal

i'm not gonna deal with python nor numpy nor any other parochial computational platform.

first, since "silence" is a perceptual property, you need to apply a weighting filter, such as A-weighting to boost the frequency components of the audio that our ears are more sensitive to and attenuate the portions we're less sensitive to.  there are other weighting curves besides A-weighting.  i like Robert Wannamaker's E-weighting the best.
then square the filtered signal (to get power or energy).
then low-pass filter the squared and filtered signal (to get an envelope).
then compare this envelope to a threshold that will correspond to the threshold of hearing.
lastly, you may need to put in a little hysteresis in the threshold comparison so that the silence/no-silence result does not jitter as your envelope crosses the threshold.

this is similar to the "gate" function in a compressor/limiter/gate process that sometimes finds expression in audio applications.  the "gate" was to make sure that something that's perceived as silence (or sufficiently low-level that it should be perceived as silence) becomes zero during the "silent" periods of time.


How to Calculate the Color Temperature / Tint of the Colors in an Image?

This paper (pdf download) gives the following formulae for calculating Correlated Color Temperature (CCT). They do not explicitly say (or I missed it), but their example leads me to infer that they are assuming RGB values in the range of 0-255.

1. Convert the RGB values to CIE tristimulus values (XYZ) as follows:

$$X = (-0.14282)(R) + (1.54924)(G) + (-0.95641)(B)$$
$$Y = (-0.32466)(R) + (1.57837)(G) + (-0.73191)(B) = Illuminance$$
$$Z = (-0.68202)(R) + (0.77073)(G) + (0.56332)(B)$$

2. Calculate the normalized chromaticity values:

$$x = X/(X+Y+Z)$$
$$y = Y/(X+Y+Z)$$

3. Compute the CCT value from:

$$CCT = 449n^3 + 3525n^2 + 6823.3n + 5520.33$$

$$\text{where }n = (x − 0.3320) / (0.1858 − y)$$

Which can be combined to form the following equation:

$$CCT = 449n^3 + 3525n^2 + 6823.3n + 5520.33$$
$$\text{where }n = ((0.23881)R+(0.25499)G+(-0.58291)B) / ((0.11109)R+(-0.85406)G+(0.52289)B)$$

I am not sure about applying this to an image, but if you just want a single, generalizing number, than you could perhaps use some sort of averaging? Either find an acceptable "average" RGB value for the image (eg. the centroid) and use that to calculate a temperature or (a much more computationally expensive option) calculate the temperature for each pixel in the image and take the average of those results.

Also, bear in mind that CCT is only an approximate metric for most colors, since only a single curve in the color space actually represents color which can be obtained from a real world black body radiator. Thus for all other colors, the calculated color temperature is simply an approximation of the black body temperature it most closely represents. Thus, for some colors (especially greens) it can actually be a somewhat meaningless value, at least in a physical sense. This is illustrated well in the following image (from the wikipedia article on color temperature).

The black line in the image represents the Planckian Locus of colors which could actually be produced by block body radiation. The smaller crossing lines represent the isotherms of the CCT approximation nearbye.

Also, since your question specifically references Adobe Lightroom, I found this while searching around:

  The sliders [in Adobe Lightroom] adjust not the black body temp of the light, but the compensation applied to the image to compensate for the black body temp of the light.  This goes the other way round.

So bear in mind that the color temperature you see on the Lightroom slider will not be the same as those calculated from the above formulae.


Identifying whether or not a cyclostationary signal is noisy or not using the cyclic autocorrelation

You're on the right track. The cyclic autocorrelation can be estimated by inverse Fourier transformation of the cyclic periodogram. If $X_T(t, f)$ is the Fourier transform of your data block (indexed by time $t$), and $\alpha$ is your desired cycle frequency, you can form the cyclic periodogram by shifting the Fourier transform to the left and to the right by $\alpha/2$ (not by $\alpha$ as you state in your question):

$$I_x^\alpha(t, f) = \frac{1}{T} X_T(t, f+\alpha/2) X_T^*(t, f-\alpha/2)$$

Then just inverse Fourier transform that cyclic periodogram. You'll have to be careful about the case where the cycle frequency $\alpha$ is not equal to an integer number of FFT bins.


Sampling: how to determine a “good” sampling rate based on a given filter

Usually, one has a noise floor requirement.  So one can set the sampling rate to over twice the frequency where the filter response (and/or spectrum in the input) goes below this desired noise floor.


Use of the harris-Moerder Nyquist Pulse Shaping Filter

Liguid-dsp implements root-Nyquist harris-Moerder (hM-3) filter using Parks-McClellan algorithm as :
void liquid_firdes_hM3(unsigned int _k, unsigned int _m, float _beta, float _dt, float * _h)


Why Is Zero padding Required for Sobel Edge Eetection?

The Sobel Filter is a $ 3 \times 3 $ matrix (it is separable, but let's ignore that).

The anchor pixel is the middle one hence to evaluate the operator on pixels on the upper row the operator needs data above them.
Same goes for pixels on the last row,  left-most column or  right-most column.

This is usually solved by padding the image (I actually use "Replication" instead of Zero Padding).

Another option is to set the output to be smaller image by not evaluating the extreme rows or columns.


Detect almost circular shapes in RGB image

The image itself is somehow blurred a little. Whole edges of trees are not visible. The only clue of tree existence seems like shadows. So I suggest using shadow information and watershed transform and local thresholding. I have done similar exercise in here.  The result is not that perfect but it will give some idea. Here are the steps. 1. Find shadows, 2. Expand shadows towards sun direction 3. Create shadow mask 4. Blur the image 5. Watershed transform 6. Threshold each watershed segment individually 7. Fit circle with Hough transform

im=imread('P0Yvl.png');
imOrig=im;

% Grayscale
im=rgb2gray(im);

% Shadows
shadow=im<60;
[xs,ys]=find(shadow);

% Grow shadow towards sun direction

% Create structural element for dilation
str_element_initial=flip(eye(25));
str_element_initial(:,1:12)=0;

shadow2=zeros(size(str_element_initial));
for i=1:7
    shadow2 = shadow2 | imresize(imrotate(str_element_initial,40-i*10),size(str_element_initial))>.1;
end

mask_shadow=imdilate(shadow,shadow2);
mask_shadow=mask_shadow-shadow;

% Plot shadow mask and initial shadow locations
figure,imagesc(mask_shadow);axis image;
hold on, plot(ys,xs,'r.')

% Blur image
im2=imfilter(im,fspecial('gaussian',[21 21],3));

% Watershed transform
L = watershed(max(im2(:))-im2);
[xw,yw]=find(L==0);

figure,imagesc(imOrig),axis image
hold on, plot(yw,xw,'r.','MarkerSize',1)

% Threshold each segment individually
tmp=zeros(size(im));

for i=1:max(L(:))
    mask=L==i;
    ind=find(mask);
    thr =multithresh(im(ind),1);
    tmp(ind)=im(ind)>thr*.8;
end

% Remove detections that don't have shadow
trees=tmp & mask_shadow;
trees=bwareaopen(trees,10);

% Find circles with Hough transform
[centers, radii, metric] = imfindcircles(trees,[5 25]);

% Result
figure,imagesc(imOrig),axis image
viscircles(centers, radii,'EnhanceVisibility',0,'LineStyle','-','LineWidth',1);


DFT independent variable as fraction of sampling rate

The DFT calculates spectral components up to $f_\mathrm{s}/2$, no matter what the input signal is. If the input signal to the DFT contains (absolutely) no spectral components between $0.33 f_\mathrm{s}$ and $0.5 f_\mathrm{s}$ the corresponding DFT outputs are zero (assuming that there is no spectral leakage).


How to obtain the probability of a detection?

You have a fairly common problem in Detection theory. Let's assume you have a $N \times N$ grid of pixel values and $N$ is small enough that the edge is a straight line in the region.  For that given line, you can formulate a hypothesis on whether and edge was on that line or not. You might also have to some other information like the amplitude of the edge and distribution of the non edge pixels.  Unfortunately, you usually don't know where the lines are, so one thing people do is something called a Generalized Likelihood Ratio GLR test.  So, if you can enumerate all the possible lines that can fall in your $N \times N$ grid and test each, you can detect it, assuming you can estimate the other parameters.  Unfortunately, that's a big search space and probably not a good approach. Even if you could do the search, since pixels will be contained in more than one search, the tests will not be independent and a closed form solution for a Detection probability will be difficult if not intractable to calculate. As an alternative, you  can do an $M$ out N test, test if there are $n$ out of $N \times N$ pixels that are part of a line.  We don't know $n$ either but if there is a single line $n$ is bound by $\sqrt{2} N$.  Again we have nonindependent tests and a difficult calculation, and the pixels that are high aren't necessarily connected as in a line.

A more practical approach is to choose $H_0$ to have no edges, and assuming we have a reasonable sampling distribution, we would look for small $p$ values. The original likelihood we mentioned would be more powerful and often people use the ideal case as an upper bound.

You might try a nonparametric test like Wilcoxon rank test on pairs of blocks.

Radon transforms are good at finding lines, so that is a possible approach as well. Picking $N$ essentially a heuristic as well.

Honestly, it all gives me a headache and I try to stick to time series.

Hopefully someone who knows better will give you a better answer


After reading your comments, the real question you're asking is:

"What's the difference between an estimator and filter?"

An estimator is a system that takes in inputs and attempts to estimate a signal (or state) given knowledge of the system. A simple example is system that takes in positional data and calculates the corresponding acceleration. In this case, your input data is position, the system is the second derivative with respect to time, and the output is acceleration. You've just estimated the acceleration by feeding it through a model.

A filter is just a system that takes an input and attempts to modify the frequency content of a signal. A Butterworth filter is a filter, not an estimator.

So the real question you have to ask your is: "What am I trying to do?"

A problem where you're trying to estimate something, you use an estimator (i.e. a system that produces an output using some physical model). A problem where you trying to change the frequency content of a signal, you use a filter. The filter doesn't care about the underlying process or signal, it's just setup to modify the frequency content of the signal.

Mathematically, a filter and estimator are the same thing (put signal in, get signal out), but they typically have different physical meanings and are used to achieve different things.

Edit: It also seems that your question may have been inspired by the 'Kalman Filter'. Even though this contains the word 'filter', its still considered more of an estimator. However, depending on how you set up the problem (i.e. how much uncertainty you have in your model and how much uncertainty you have in your feedback), the system will experience different frequency response characteristics.


How to reconstruct the signal using overlapping frames in MATLAB

unbuffer in Python; helpful post. Nice visuals in MATLAB docs (bottom).


First, I suggest you take a look at the STFT.

Let say you use $N$-point FFT, then you filterbank has $N/2$ filters (assuming real signal). Usually, you divide your signal into $K$ segments, each $N$ samples long (could be overlapping) and apply the FFT to each segment. Every time you apply the FFT you get a single sample for each filter. Therefore your filterbank has $N/2$ filters with $K$ samples.


my experience in audio is that for IIR filters, there is usually a cascaded Second-order Section (SOS) for every little feature (a bump or an edge in the frequency response deliberately placed there by the user).  often just a single second-order IIR is used for many little jobs where some frequency is given a little boost or cut.  we call those IIR filters a Parametric Equalizer.

however, for a fixed and sharp feature (like a brickwall of some sort), then the order of the IIR filter will be related to the sharpeness of that feature.  for a brickwall, there is a 6 dB per octave of steepness in the wall for every integer order of the filter.

for FIR, Fat's answer is good.  the order of the filter is 1 less than the number of taps of the FIR.  the minimum order again is related to the sharpness of some particular feature in the frequency response.  there are a couple of heuristic formula for initial guessing at what order you will need.  i will try to express them in a way to compare apples-to-apples.

using a Kaiser-windowed FIR design the filter order (one less than the number of taps) is about:

$$ N \ge \frac{(A-7.95)F_\text{s}}{14.357 \cdot \Delta f} $$

where $A$ is the amount of attenuation in dB of the brickwall, $\Delta f$ is the frequency width of the transition between the nominally un-attenuated and attenuated portions of the frequency response, and $F_\text{s}$ is the sample rate.

using Parks-McClellan will reduce $N$ a bit:

$$ N \ge \frac{(\tfrac12 A + \frac{4.343}{r_\text{p}} - 29.7)F_\text{s}}{14.602 \cdot \Delta f} $$

where the new symbol $r_\text{p}$ is the passband ripple in dB.

this appears to be twice as good as i was expecting.  so i am not sure i did this right, but i was being very careful with the equations from O&S about the P-McC alg.


Homography valid for patch within the image but not the entire image?

Are you sure that the points are all in front of the camera (z > 0), if you apply a a projection like $(x, y) = (x/z, y/z)$ then points in a 3d quadrilateral where vertices 0 and 1 have negative z, you could end with a 2d figure like you described


Having converted a grayscale image to binary using `im2bw()`, is the converse operation possible?

As you've stated the problem, you can't. Converting a grey-scale image to black and white (binary) you're comparing pixel values (0-256 for an 8-bit image) to a threshold (that you define) and giving each pixel a value of 1 if the pixel is above the threshold, and 0 if it's below the threshold. The binary image doesn't contain information about by how much the pixel value was above or below the threshold.
If you want to get the pixel values from the original image which are above the threshold, you can multiply the binary image and the original image. All the values in the original image below the threshold will be zeroed out when they are multiplied by the zeros in the binary image. All the pixel values above the threshold in the original value will be multiplied by one and will retain their original values.
Unfortunately, if you're given an arbitrary binary image with no information about the threshold used, or the original image, you cannot recover the original grey-scale values.
Hope this helps.


Why is it called 'multi-level thresholding', if only the best threshold is selected?

Using k1 and k2 you can calculate a between class variance value using the Otsu formula. However, there may be many combinations of k1 and k2 for which you have to calculate between class variance values to see which of them gives the highest between class variance value. That's what 'max' means. So whichever k1 and k2 gives you the highest between class variance, you choose those values of k1 and k2 as your thresholds. Same can be done if you want more thresholds.


Demodulation of a BFSK frequency hopping signal using FFT

At some point, before the fft() function is called, the sampling frequency must be known (from ADC? or other resampling functions?)... use that to find the frequencies of the bins as you would in Matlab.


System Identification using LMS Adaptive Filter

Because an LMS estimator will, over time, "average out" uncorrelated zero-mean noise. It's pretty much in the name.

But yes, you're right, there is a noise component in the estimate; one of the qualities of an estimator is how little the noise variance influences the estimate variance after a given length of observation.

That's the case, however, for all estimators: you measure signal + noise, and you estimate parameters from that. The parameters must be somewhat noisy; otherwise, there's something broken with your noise model.


What is the importance of signal bandwidth in radar systems for velocity and localization estimation

Bandwidth affects the "gate/cell width" of your measurements. If you think about an LFM waveform reflecting back off an object and going through a matched filter to determine the range, a signal with a very low bandwidth is going to have a more ambiguous return time (range) than a more aggressive chirp. The rough calculation for the smallest range resolution for a certain bandwidth is $\frac{c}{2B}$.
I imagine the relationship is reversed for doppler processing. A signal with low bandwidth would provide a lot of doppler resolution because there isn't a great deal of ambiguity as to the return's center frequency, so the doppler gate width should increase proportional to bandwidth. That is a guess, though.
Radar and other signal processing disciplines involve these tradeoffs a lot. Certainty in frequency versus certainty in time and vice versa, maximum  unambiguous range versus improvement factor from pulse repetition frequency, etc.


Using a Wiener Filter to Estimate a Transfer Function

The desired MSE is application dependent, so there can be no general
rule. If the approximation doesn't satisfy your needs you can
increase the filter length to obtain a better match.
There is no straightforward way to change the FIR Wiener filter solution to an IIR solution because the IIR formulation results in a set of nonlinear equations which have no closed-form solution. The IIR solution might also be unstable, so FIR filters are a much more practical choice when computing a Wiener filter.