Otsu's method problem? I have a problem with Otsu's method. I would like to implement histogram to compute how many ones' and zeros'. For Otsu's method, I have a code from a book and it is running without problem. This is the picture

.

I've implemented a code to compute the histogram. This is the output of my code.

------------------
    0 : 14332
    1 : 2102
    2 : 1387
    3 : 745
    4 : 221
    5 : 117
    6 : 27
    7 : 9
  248 : 6
  249 : 12
  250 : 48
  251 : 162
  252 : 396
  253 : 830
  254 : 1377
  255 : 28554
------------------

Here's my code:

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

This is my C++ code.

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

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

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

        }
    }

    return histogram;
}

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

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

Why this is the result? What I'm expecting is how many 255 and 0 (black and white only). Moreover, if I run another picture, I get another results with the range 0-255.

jpg is a compressed storage format and no matter how binary your image is, you will almost always end up with gray values in the resulting image (due to compression).
Please save it as png and try again.

(Also note that you still have 0s and 255s as the dominant bins, because jpg generally degrades around the edges, which constitute a small portion in the image.)


Serial to Parallel Converter in many texts about modulators I have seen there is a component called "Serial to Parallel Converter". For instance, let's consider this scheme and let's focus on the second block (reference):

The text about that block says:

  The input serial data stream is formatted into the word size required for transmission, e.g. 2 bits/word for QPSK, and shifted into a parallel format. The data is then transmitted in parallel by assigning each data word to one carrier in the transmission.

Somewhere I have read also that this device converts a signal from serial to parallel by halving the transmission rate.

But,

I do not understand what serial to parallel conversion exactly means, and if I have to see it from and electronic and circuital point of view (from a signal between a wire and GND we get two signals between two wires and GND) or from a signal theory point of view (but I do not know how to see it).
Which is the math relationship of the output signals (between them and between them and the input signal)?

So, I myself have taught based on material that uses that scheme, "P/S" and "S/P" after and before the transforms.

Personally, it's nonsense.

What the author tries to say is:

  The IFFT is a mapping of sample vectors to vectors. So, you need a vector as input, not a stream of samples.

What they instead say is:

  We use the terminology from very basic digital logic to deal with complex values.

That's not only a bit awkward, if you asked me, but also inaccurate:

modern FFT implementations in hardware actually do take in stream data, so that operation isn't there in practical hardware implementations
in software implementations, you basically never even deal with samples coming in serially – they always appear en block, in some memory location, so nothing's "serial" to even begin with
S/P nor P/S even say that you should be fully accumulating one vector of $N$ samples, then move on to the next vector, and have zero overlap between these – when applied to FIFOs (which is where the term comes from), that's usually not how they work.

So, simply think of S/P as "get $N$ samples, present these $N$ to the next block as vector".

And P/S is "take this $N$ long vector, and give one sample after the other".

So, these blocks do exactly nothing to your signal - it's just a reinterpretation, if you will, between things that are logically "one after the other", and things that are logically "a vector of $N$ values".


Can you decimate / downsample a signal in frequency domain just like you can interpolate / upsample it? To interpolate a signal I can just zero pad it in the frequency domain. If I want to decimate the signal, can I just discard some part of the frequency domain?

So in MATLAB this works:

sInterpolated = abs(ifft(ifftshift(padarray(fftshift(fft(s)), nextra/2, 0, 'both'))));

Would the opposite also work (correctly)?

Sshifted = fftshift(fft(s));
Sdecimated = Sshifted;
Sdecimated(1:nless/2) = [];
Sdecimated(end-nless/2+1:end) = [];
sDecimated = abs(ifft(ifftshift(Sdecimated)));

Here s is the (real) input signal, nextra is the number of extra samples in the interpolated signal and nless is the number of decreased samples in the decimated signal.

Your interpolation operation is an upsampling operation with subsequent low-pass filtering. This translates to zero-padding in the frequency domain. Given your original signal $x[n]$, you upsample it to become

$$y[n] = \begin{cases}x[n/M]& n/M \in \mathcal{Z} \\ 0 & \text{else}\end{cases}$$

to become an M-times upsampled signal. Then, your interpolated signal is

$i[n] = h[n] * y[n]$ where $h[n]$ is a Dirichlet kernel/periodic sinc function with zeros of distance $M$. In the spectral domain, this corresponds to

$$Y(f) = \sum_{m=0}^{M-1}X(M(f-m))=X(Mf)$$
I.e. a M-times periodic repetition of its spectrum (Note that a periodic spectrum creates the discrte time-domain signal).  Here, I assume $f=0...1$ for the discrete case and the last equality follows from periodicity of $X(f)$ here. Also note that your original spectrum is narrowed, since your sampling frequency increased.

Now, you filter out the spectrum repetitions by multiplying with a rect-function in the frequency domain. This corresponds to convolution in time with the sinc function.

Now, to do the reverse (i.e. decimation), this is what happens: Let $x[n]$ be your original signal you want to decimate. You calculate

$y[n] = x[Mn]$ and this corresponds to the following in the spectrum:

$$Y(f) = \sum_{m=0}^{M-1} X((f-m)/M)$$

I.e. your original spectrum is repeated M times overlapping each other and summed together. This in general creates aliasing. Hence, a simple removal of higher frequencies in the DFT domain is not the same as taking every Mth sample from the time domain signal. Instead, you need to take the spectrum, sum it up repeatedly (such that it becomes M periods within your full frequency band) and finally perform the ifft of one period of the spectrum.

Your approach works, if your signal is bandlimited enough, if it does not create aliasing, i.e. the bandwidth of the signal is smaller than 1/M. Then, the shifted spectra dont overlap and you can just remove do ifft of the smaller portion. Alternatively, you can see your approach as filtering with an aliasing filter (i.e. lowpass with cutoff 1/M) before downsampling.


strange artifact in depthmap of a mirroring glass I hope this the right place for this question, not sure.
I have a Asus Xtion Pro live camera, so an IR structured light sensor, on a shop window. When the shop window glass mirrors (there is more light inside than outside the shop), I see this mirrored image also in the depth map. It's clearly an artifact, but how is it possible?

Ignore the holes in the depthmap (black areas) for the moment, it's the not the focus of this question, but look at the red circle, this is a table, that is located inside the store (as you can see in the rgb frame), it should not be present in the depthmap.
As far as I know a structured light sensor project a light pattern on the scene and through a camera read how the pattern is deformed on the different surfaces (near, far) to calculate the depth (distance of every point in the detphmap). Of course the table doesn't exist in the scene, so should not exist in the depthmap too! Any idea?

Structured or not, light gets reflected and the sensor has no way of telling if the scene is imaged directly or by reflection.

To try and reduce the effect of reflections from the glass, you might want to use a polariser. Photographic polarising filters have lots of useful applications and one of them is reducing glare when you are trying to shoot artifacts within a museum glass case (for example).

You can find really cheap polarising filters in the second-hand boxes in camera shops, or buy polarising sheets which you can cut into various shapes.

One thing that would be good to note however is that the camera may already be using a polarising element internally which might interfere with an additional polariser which will worsen the situation.

Hope this helps.


Sequential Non Linear Least Squares Problem I have the the following non-linear function,
$$s(x;A_k,\mu_k,\sigma_k)=\sum_{k=1}^2 A_k \exp\left(\frac{-(x-\mu_k)^2}{\sigma_k^2}\right)$$
with unknown (but deterministic) parameters $A_k,\mu_k,\sigma_k$.

The received signal is as follows: $$r(x;A_k,\mu_k,\sigma_k)=t\cdot s(x;A_k,\mu_k,\sigma_k)+n(x)$$ where $n$ is standard white Gaussian noise and $t\geq 0$.

The task is to estimate/update the model parameters $A_k,\mu_k,\sigma_k$ online, when the next data point arrives, because the data points arrive sequentially.

Here the data points DO NOT ARRIVE randomly from the model. They arrive in the order as given by the model. Namely, the incoming points are first around zero level, then they increase, make a peak, and then decrease, and the same happens for the second exponent.

We know that $\sigma_1^2\approx \sigma_2^2$ and $\mu_2-\mu_1$ allows that two distinct peaks can be visible. The incoming values are allowed to be evaluated over a moving window of size $5\%$ or less of the whole data. All data points lying behind this window are simply lost, due to memory constraints.

I was thinking about using an iterative non-linear least squares algorithm but this method doesn't consider about how the data points arrive. It assumes that data points arrive randomly. However, in my case, with a few data points from the first exponent, if I know that those points are from the left tail, for example, and the noise level is low, I can even almost completly know $A_1$, $\mu_1$ and $\sigma_1^2$ but this is never possible with the regular non-linear least squares. Hence, my question is the following:

  What is the best method to solve this problem?

I have the the following non-linear function [...]

This would be a bimodal distribution with two different means and deviations ($\mu_1, \sigma_1, \mu_2, \sigma_2$), two sub population sizes $A_1, A_2$ and an extra parameter $t$ that modulates $A_1, A_2$.

  Here the data points DO NOT ARRIVE randomly from the model. They arrive in the order as given by the model. Namely, the incoming points are first around zero level, then they increase, make a peak, and then decrease, and the same happens for the second exponent.

What you are describing is plain simple curve fitting of some curve to data points as opposed to fitting data to a distribution.

The typical fitting scenario involves the use of Newton's method to find those parameters that minimise the value of the error function and here it is possible to differentiate the sum of the Gaussians with respect to the parameters (which is a common requirement for both linear and non-linear methods).

Here are a few practical examples: One, two, three (different platform).

If you can use details of your problem to fix the starting points of the fit, it will help immensely with convergence. There are some cases for example where the method converges but when you look at the fit, one curve has simply absorbed the other.

Can it be done with an online algorithm? No. Not in this case.

In general, there are online algorithms that will "learn" the coefficients of a model in a "sample-by-sample" basis. But this refers to working out parameters of a model that minimise an error function. Evaluating the error, implies a full "snapshot" of what is being fitted (i.e. all of the data to obtain one estimate of $s$).

Here, it is impossible to predict what sample comes next and from which possible realisation of the 6 parameter model. But there might be a thin region of cases where you could make an inference provided that you can apply additional constraints.

You might see this problem as a combination of a tracking and a classification problem where you can say "I know how the curve is supposed to evolve and all I have to do is track its evolution and raise a flag when the uncertainty falls below a level", but then again you don't really know which 6 parameter realisation to track. This is what you are trying to infer.

It is impossible to "guess the future" unless there are certain constraints.

What do the first few "... incoming points are first around zero level,..." tell you about the value of $A_1,\mu_1, \sigma_1, A_2, \mu_2, \sigma_2$? Nothing. You can still attempt a "fit" (batch or online) but you will get back nonsensical results and a huge error. So, we have to wait for more data. But then again, if we wait long enough to be sure we have one of the curves "in", then, we probably can make a quick guess about the $A$s provided that we know what $A_1+A_2$ is in general (another constraint). But still, this leaves us with no way of guessing what the $\mu_2, \sigma_2$ would be. We still have to wait for more data.

If the problem you are dealing with can already inform you about nonsensical values for the 6 parameters or at least the range they are supposed to be in or some sort of relationship between the two triplets of parameters, or anything else that could be used to reduce your search space, then, you might be able to pre-calculate a very large number of realisations of $s$ for all the different possible 6 parameter combinations and then track the evolution of each one as a new sample comes in. Which is basically, nearest neighbour classification.

But still, you would probably need to wait almost for the whole curve to come "in" to be able to make a certain decision.

If you start tracking a single Gaussian curve with your eye, you cannot really tell what $A, \mu, \sigma$ are until the trace starts curving to begin forming the peak. Add to this the noise and you have more than one possible curves that could be passing from the same data. So, you probably have to wait for at least $\frac{3}{4}$ of the data to come "in", before you can make a valid estimate.

Hope this helps.


DFE for QPSK in MATLAB I need to do a DFE for QPSK simulation in MATLAB.
The system consists of a QPSK signal transmitted with power = 1 which is then pulse shaped with a square raised-cosine filter. Add AWGN and a filter receptor with ISI = [0.25 1 0.25].

Summarizing and clarifying, the sequence is:

QPSK signal -> Upsample -> Rcosine -> AWGN -> ISI [0.25 1 0.25] -> Downsample -> DFE -> BER counting.

I have coded the algorithm, with guidance from Dig. Comm. Lee-M., but I am thinking about problems with the Feed-Forward filter. The objetive is to obtain a signal with NO ISI, and the correct BER for a given SNR.

The code is:

% DFE.
tapsffe=7;
tapsdfe=8;
% Gain
betaf=0.0001;
betab=0.0001;
% Coefficients init. Ck for FFE and Dk for DFE
ck=zeros(1,tapsffe);
ck(7)=1;               % Last tap for FFE=1
dk=zeros(1,tapsdfe);

% Signal to filtering
rk=zeros(1,tapsffe);  % Signal to filter in FFE
ak=zeros(1,tapsdfe);  % Signal to filter in DFE

% Filtered signal, and correspondings error
n=length(s6);
s7=zeros(1,n);
e7=zeros(1,n);

for k=1:n

    rk = [s6(k), rk(1:tapsffe-1)]; % FILTERING in FFE

    % Output from Feed-Forward and FeedBack
    ykf = sum(rk.*ck);
    ykd = sum(ak.*dk);
    yk = ykf-ykd;

    % Slicer for QPSK. sqrt(2) for POWER = 1.
    dec = (-1*(real(yk)<0)+1*(real(yk)>=0) + 1i*(-1*(imag(yk)<0)+1*(imag(yk)>=0)))/sqrt(2);

    % Obtained value for filter ecualized signal
    s7(k) = dec;

    % Error:
    error = dec-yk;
    e7(k) = error;

    % Coefficient adaptation
    ck = ck - betaf*error*conj(rk);
    dk = dk + betab*error*conj(ak);

    % FILTERING for Feedback filter with SLICER decisions
    ak = [dec, ak(1:tapsdfe-1)];
end

Apparently I have problems with the FFE, because I do not see ck taps to accommodate to the inverse of the channel response with ISI. I have little doubt of its coefficients adaptation. Signal for BETA, and adaptation for ck with the corresponding rk, or the fliplr coefficient...

If we focus on finite-duration impulse response filters and apply the
MSE criterion to optimize the coefficients of the feedback and feedforward filters.

The equalizer output can be expressed as:
\begin{equation}\hat{I}_{k}=\sum_{j=-K_{1}}^{0} c_{j} v_{k-j}+\sum_{j=1}^{K_{2}} c_{j} \tilde{I}_{k-j}\end{equation}

where $\hat{I}_{k}$ is an estimate of the $k$ th information symbol, $\left\{c_{j}\right\}$ are the tap coefficients of the filter, and $\left\{\tilde{I}_{k-1}, \ldots, \tilde{I}_{k-K_{2}}\right\}$ are previously detected symbols. The equalizer is assumed to have ( $K_{1}+1$ ) taps in its feedforward section and $K_{2}$ in its feedback section.

The MSE criterion result in a mathematically tractable optimization of the equalizer coefficients; based on the assumption that previously detected symbols in the feedback filter are correct, the minimization of MSE:
\begin{equation}J\left(K_{1}, K_{2}\right)=E\left|I_{k}-\hat{I}_{k}\right|^{2}\end{equation}
This leads to the following set of linear equations for the coefficients of the feedforward filter:

\begin{equation}\sum_{j=-K_{1}}^{0} \psi_{l j} c_{j}=f_{-l}^{*}, \quad l=-K_{1}, \ldots,-1,0\end{equation}
where
\begin{equation}\psi_{l j}=\sum_{m=0}^{-l} f_{m}^{*} f_{m+l-j}+N_{0} \delta_{l j}, \quad l, j=-K_{1}, \ldots,-1,0\end{equation}

The coefficients of the feedback filter of the equalizer are given in terms of the coefficients of the feedforward section by the following expression:

\begin{equation}c_{k}=-\sum_{j=-K_{1}}^{0} c_{j} f_{k-j}, \quad k=1,2, \ldots, K_{2}\end{equation}

The values of the feedback coefficients result in complete elimination of intersymbol
interference from previously detected symbols, provided that previous decisions are
correct and that $K_2 > L$


Should I use CV_HAAR_SCALE_IMAGE while using LBP CascadeClassifier? I trained a cascade-LBP to detect Lollipops with 1000 images, now I'm trying to "adapt" the openCV HAAR-Cascade example to use my LBP trained .xml but I'm not sure about the "HAAR_SCALE_IMAGE" flag, I read that optimizes the downscale to increase the HAAR features is it right to use it with LBP?

//-- Detect Lollipops
  lollipops_cascade.detectMultiScale( frame_gray, lollipops, 1.1, 2, 0|CV_HAAR_SCALE_IMAGE, Size(30, 30) );

An example of imput image:

For more info:
http://docs.opencv.org/2.4/doc/tutorials/objdetect/cascade_classifier/cascade_classifier.html

The OpenCV documentation suggests using CASCADE_SCALE_IMAGE when performing a detectMultiScale with LBP features.

lollipops_cascade.detectMultiScale(frame_gray, lollipops,1.1, 2, 0|CASCADE_SCALE_IMAGE, Size(30, 30));

Even though, it works perfectly with  CV_HAAR_SCALE_IMAGE flag instead of CASCADE_SCALE_IMAGE.


Suggested Preprocessing methods for OCR on Circular Images Hello this is my sample image I am going to do real time character detection on images like that.

I've tried SURF, SIFT, MSER and template matching on original image without any preprocessing. I can detect characters but when the image changes this feature detection methods does not work because there isn't a lot of interested points in the characters on the image. I am not exactly sure what kind of preprocessing I should do for higher accuracy.

so my question is what methods can you guys suggest me to do for character detection with higher accuracy.

thanks in advance.

Here is what I did for a client (What you are asking is the same).

Assuming that you have access to certain type of a pattern on the image (or the center of the hole), you could always detect the template to obtain the location of a possible unwarp:

Note that in the transformed image, two region of interests are defined and the region within which we would like to read the characters is the torus lying between two circles. Using this information, we could warp this image using an inverse polar transform to obtain:

I also provide the MATLAB code for this transformation at this file-exchange post.

In the next step, thanks to the template matching, the characters to be read always lie in the local region defined in the image above. Now, if you have appropriate illumination the character segmentation is nothing but a trivial set of image processing operations:

The connected components obtained eventually are fed into a neural network or an SVM, first to train and then to classify. This part is composed only of a basic machine learning task, which I won't go into the details, here.

Don't forget to take good care of the illumination, as it is almost solving the problem half way.


Time variant and Time Invariant Systems I was trying to solve the below system is a Time variant or Time Invariant System:

$y(n)=nx^{2}(n)$

Now for a system to be time invariant the system property should be :

y(n,k) = T[x(n-k)]

OR

y(n,k) = y(n-k)

for all possible values of k.

Coming back to system equation:

$y(n,k) = T[x(n-k)] = nx^{2}(n-k)$ delayed by k sample

And

$y(n-k) = nx^{2}(n-k)$

Hence

$y(n,k) = y(n-k)$

The the system is time invariant, but the solution in book states that the system is time variant. I am not sure where I am going wrong or my understanding is not correct?

Here is a simple way of looking at it which avoids any algebra and uses minimal symbolism. Let us calculate the output $y$ for an input $x$ which is a unit pulse function. We have that

$y[-2] = (-2)(x[-2])^2 = 0$ since $x[-2] = 0$.
$y[-1] = (-1)(x[-1])^2 = 0$ since $x[-1] = 0$.
$y[0] = (0)(x[0])^2 = 0$ (regardless of the value of $x[0]$).
$y[1] = (1)(x[1])^2 = 0$ since $x[1] = 0$.
$y[2] = (2)(x[2])^2 = 0$ since $x[2] = 0$.

and so on.  The entire input is blocked and the output $y$ is always $0$.

What about the output $\hat{y}$ when the input $\hat{x}$ is just the previous $x$
delayed by one time unit? The only difference now is that $\hat{x}[1] = 1$ and
we have that

$\hat{y}[1] = (1)(\hat{x}[1])^2 = 1$ since $\hat{x}[1] = 1$.

It is, I hope, obvious that $\hat{x}[n] = 0$ for all integers $n$ except $n=1$.
So, since the input $\hat{x}$ is just a delayed version of input $x$, is
is true that the output $\hat{y}$ is just a delayed version of output $y$?
What, if anything, does this tell you about the system properties?


Filtering 50Hz using a notch filter in MATLAB I've read up a lot about this, but haven't been able to piece everything together successfully, so I'm looking for some help.

I need to filter 50 Hz from a signal. It looks like the best options are either a notch filter or a LMS filter, but I don't have a copy of the noise so a notch filter seems to be the best choice.

I don't need the frequencies normalised, as I know the sampling frequency (16kHz), and the duration is 30 seconds. The bandwidth can be fairly tight, 49.5Hz ~ 50.5Hz should be fine.

It looks like I need to use a combination of filter and iirnotch, but I'm not entirely sure how.

If someone can bring this all together I would greatly appreciate it. Thanks.

I'm not sure what iirnotch does, but this is how to design the notch filter by hand.

fs = 20000;             % sampling rate
f0 = 50;                % notch frequency
fn = fs/2;              % Nyquist frequency
freqRatio = f0/fn;      % ratio of notch freq. to Nyquist freq.

notchWidth = 0.1;       % width of the notch

% Compute zeros
notchZeros = [exp( sqrt(-1)*pi*freqRatio ), exp( -sqrt(-1)*pi*freqRatio )];

% Compute poles
notchPoles = (1-notchWidth) * notchZeros;

figure;
zplane(notchZeros.', notchPoles.');

b = poly( notchZeros ); %  Get moving average filter coefficients
a = poly( notchPoles ); %  Get autoregressive filter coefficients

figure;
freqz(b,a,32000,fs)

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


Why FIR uses often direct-form instead of transposed-form? Even though tranposed-form has less computational burden and memory requirements, why FIR uses mostly direct-form?

I don't see how memory or computational burden can possibly be reduced between the two implementations (note this is not implying simplicity, just number of computations required, please read on...) Usually the decision to use transposed-form is for high speed FPGA implementations with a large number of taps as you can eliminate a long adder tree which would otherwise limit the maximum clock frequency to a lower number due to the accumulated effect of multiplier and adder delay, but this can get more complicated for a fixed precision implementation since the bit width must grow as you proceed in the filter. (This consideration is particularly important in IIR implementations where feedback exists but affects FIR implementations with large adder trees, and also a reason to factor longer IIR filters in cascade form). At lower speeds the direct-form is more convenient to implement as you can use a single extended precision accumulator.

I include a summary of the basic forms applicable to FIR/IIR below for benefit of readers that may be less familiar, and see link from Ben in the comments for more advanced structures specific to FPGA implementations:


Convolution theorem for cross-correlation Forgive me is this is an ill-posed question.

Is there any such thing as a 'convolution theorem' for the cross-correlation.

Namely, the convolution theorem states that:
$$
x[n] * h[n] = \sum_{k=0}^{\infty}h[k] x[n-k] = X(f)H(f)
$$
Where $X(f)$ is the Fourier transform of $x(t)$.

Is there any equivalent theorem for cross-correlation, i.e.:
$$
\textrm{corr}(x[n],h[n])=r_{xh}[k]= \sum_{k=0}^{\infty}h[k] x[n+k] = F_?[X(f),H(f)]
$$

In other words, is there a way to express $r_{xh}[k]$ in terms of $X(f)$ and $H(f)$, or even in terms of the autocorrelation $r_x[k]$ and $r_h[k]$.

I tried self-deriving this but couldn't come up with anything.

Thanks!

To stick to proper notations, if $\mathcal F \left\{x[n]\right\}=X(\omega)$ and $\mathcal F \left\{h[n]\right\}=H(\omega)$, what you mean is the Fourier transform of a convolution gives the product of the transforms of each term:

$$
\mathcal F \left\{x[n] * h[n]\right\} = \mathcal F\left\{\sum_{k=0}^{\infty}h[k] x[n-k]\right\} =X(\omega)\cdot H(\omega)
$$

As already mentioned, there is a similar formulation for correlation, the correlation theorem. And again, the Fourier transform of the correlation:
$$
\mathcal F\left\{\mathrm{corr}\left(x[n],h[n]\right)\right\}=\mathcal F\left\{r_{xh}[n]\right\}= \mathcal F\left\{\sum_{k=0}^{\infty}h[k] x[n+k] \right\}= X(\omega)\cdot H^*(\omega)
$$

The product $X(\omega) H^*(\omega)$ is called cross-energy density  spectrum and noted $S_{xh}(\omega)$, and this is just $\mathcal F\left\{r_{xh}[n]\right\}$. The $^*$ sign is for conjugate.

Note that your convolution/correlation are on discrete-time signals $x[n]$ and $h[n]$. So $X(\omega)$ and $H(\omega)$ are transforms of $x[n]$ and $h[n]$ respectively, and not of continuous-time signals $x(t)$ and $h(t)$.


Downsampling and low pass filtering in one step? Image pyramids are created by applying a lowpass filter (Gaussian) and then decimating the image (keeping only every n'th sample).
In The Pyramid as a Structure for Efficient Computation the start image is 257x257. Then the next level is 129x129, 65x65 etc.
Without decimation, one could get the same image size by applying a convolution of stride 2 and padding 2. For $i = 257$, $k = 5$, $s = 2$ and $p = 2$ the result is $\lfloor\frac{i + 2p - k}{s}\rfloor + 1 = 129$. The kernel is still Gaussian so a low pass filter.
In Convolutional Neural Networks it is actually quite common to use a higher stride to reduce the image size. This is an alternative to max pooling where one decimates the image by computing the maximal value in $k\times k$ image patches (obtained by a filter).
My question is now: Is it possible to combine decimation and low pass filtering in one step? Not necessarily only for images but also for general signals.

Is it possible to combine decimation and low pass filtering in one step? Not necessarily only for images but also for general signals.

Yes, that's what people usually do when they implement downsampling: since of the output of the anti-aliasing filter, you throw away N-1 samples, why even calculate these?
The trick is to decompose your filter into polyphase components, which enables you to run the resulting filter operation only once per output of the downsampling, instead of once per input. There's plenty of reference implementations - from GNU Radio's decimating FIR filters, to rescalers in image processing hardware.
Think of it this way:
The trick is to take your original filter $[h_0, h_1, h_2, h_3, \ldots, h_N, h_{N+1}, h_{N+2},\ldots,h_{2N}, h_{2N+1}, \ldots]$ and just split it up into filters where there's only one non-zero entry every $N$ coefficients. Choose the non-zero-value positions so that the first polyphase component filter gets $h_0, h_N, h_{2N},\ldots $, the second gets $h_1, h_{N+1}, h_{2N+1},\ldots$ and so on.
Add up the result of these filters, when you feed in the same input, to "undo" the splitting. This doesn't change anything, it's the same filter, just split into $N$ filters with many zeros in them, but with the non-zero elements in different positions.
After the addition, you decimate by $N$. Ok, you can do that before the addition, so now you have one input stream, fed into $N$ subfilters, each with a lot of zeros in them, each followed by a decimation by $N$.
Now you have a special kind of filter that only had every Nth filter tap occupied, so the first subfilter's coefficient vector is $[h_0, 0, \ldots, 0, h_N, 0, \ldots, 0, h_{2N}, 0 \ldots]$, and you'd decimate by $N$ afterwards, you could as well just swap decimation and filter, and just use the filter $[h_0,h_N,h_{2N},\ldots]$. The two things are identical in effect; this is called Noble Identity.
So, we can "pull the decimation up front" for that filter. You can, in fact, do that to all subfilters (you will have to add delay so that it works out mathematically for the non-zero-phase polyphase components, but the idea doesn't change. You have one input stream, going into $N$ different delays, decimate-by-$N$ decimators, subfilters, and a summation.
As it happens, this means that only one "branch" at a time actually gets input per input cycle.


Data Adaptive technique I am doing data analysis. I used the wavelet transform and now I am trying the Hilbert–Huang transform (HHT). In the literature, I read that Hilbert–Huang transform (HHT) is an adaptive technique. I tried to find how it is adaptive or what is the meaning of adaptive in this context. I didn't get any answer.

The Hilbert-Huang signal decomposition possesses a non-adaptive part (the Hilbert transform) and an adaptive, or data-driven part: the decomposition of  a signal as a sum for Intrinsic Mode Functions (IMF), that are obtained by a deflation process, using extrema and envelopes of the signal, recursively.

Thus, IMFs are adaptive with respect to the signal: if you change the value of one sample, the IMFs can be drastically different. For some classes of processes, namely fractional Brownian motions, the
Empirical mode decompositions [can act as] as data-driven wavelet-like expansions. In other words, under some conditions, an adaptive algorithm may behave like a non-adaptive one.

There is one code for Marginal Hilbert Spectrum on Matlab Central.


Digital Gaussian White noise signal generation in C++ I am trying to generate a band limited noise signal on a Raspberry Pi. I calculate samples at every 1/48000 sec on the Pi. If I generate normally distributed random numbers as data points at every sample would it be gaussian noise signal. Here is what i am having a hard time to understand, how to interpret the frequency of the noise signal. And also how to limit the frequency values to be between a band?

If you are generating sample of white noise (Normal distribution value) with frequency of 48 kHz, then you already have this in frequency range of 0 - 24 kHz (at least you should if generator is truly random). If you use same samples with different sampling frequency, then they will also be a white noise.

Now if you want to make it a band-limited white noise, then you must use a digital filter. Choice is up to you, but for sake of simplicity I would recommend FIR filter. You must use it to filter in real time your incoming samples. How to do it? Well, that's totally up to you, you can even use short buffer to perform filtering by linear convolution for each generated sample but it won't be very efficient.


Real-Time peak detection (amplitude) I capture signal (sound) in 44100hz and I need get response if there is peaks at amplitude, actually. Is this possible? How I can do this in real-time?

On a modern Pentium computer, it is possible to do many algorithms in real-time at 44.1 kHz. I suspect that what you want to do is possible, but what isn't possible is to tell you how to do it without a better explanation of what you are trying to do, which should include a better definition of what you mean by "peaks at amplitude".


When to Apply Circular Convolution Formulas? Context

I am studying the family of Discrete Trignometric Transforms (DTT): Discrete Cosine Transforms (DCT) and Discrete Sine Transforms (DST). And trying to understanding more their properties, I came up to the following concepts.

Linear Convolution

$$y[n] = x[n] \star h[n] = \sum_{k=0}^{N-1}x[k] \, h[n-k]$$

Finite Sequence

$$x = {x_0, x_1, ..., x_{n-1}}$$

Which is usually extended by:

$$
x[n]=
\begin{cases}
x_n, &\text{for } n \in [0, N-1]\\
0, &\text{o.w.}\\
\end{cases}
$$

Then, we have this section of this book:

Question

  How to properly use the formulas for circular convolution and skew-circular convolution? Why to use these if they are equivalent to the traditional convolution formula?

My first thought was for avoiding re-calculations. But I didn't quite visualized and example.

[Updated]
I know that, as defined by the convolution theorem property, the convolution in frequency domain, (for instance the Fourier Domain) requires less operations for signals above a certain size, then operating it in the spatial domain. This is not what I am looking for.

Circular convolution can be done using FFTs, which is a O(NLogN) algorithm, instead of the more transparent O(N^2) linear convolution algorithms.  So the application of circular convolution can be a lot faster for some uses.

However, with a tiny amount of post processing, a sufficiently zero-padded circular convolution can produce the same result as a linear convolution, while still being a lot faster when using FFTs.  This is because the tail portion of a sufficiently long zero-padded convolutional result is all zeros, rather than being a non-zero tail result that mixes/sums with the beginning of the convolution result when doing circular convolution, which can make a mess out of non-perfectly-periodic inputs.

For sequences of windows of data, one can extend this to overlap-add or overlap-save FFT fast linear convolution.


How can I graph the amplitude for a given frequency window over time? I have simple audio files of single plucked notes from a guitar:

I want to graph the amplitude of each partials over time, rather than just have this color-coded "dB" intensity graphing I have here.

I can delete the extraneous frequencies in Izotope RX easily enough which I think does the same thing as FFT windowing:

Is this the right approach? How can I next get the amplitude of this window over time?

I know basic C++/JUCE but mostly for synth design. I'm not a pro on DSP, FFT, etc.

I appreciate any guidance.

Do you have access to MATLAB or Octave or similar?

In general, I would:
 - Compute the spectrogram (STFT) of the signal you're interested in.
 - Perform any required scaling to get a good amplitude estimate.
 - Apply some thresholding to treat the noise floor.
 - Devise a scheme for identifying the amplitude maxima of the harmonics you're interested in.

This final point will be much easier if you know the fundamental. For example, if you're playing a note at 440Hz (A440), and you know you only want to examine the fundamental, second, and third harmonics, you can easily define some brackets to search within.

If you don't know the fundamental, you will have to implement some method for determining what the fundamental is. However, assuming you do know the fundamental, I've provided a MATLAB example.

This example first builds an arbitrary signal with some sinusoids that decay at differing rates, then analyses the signal to identify harmonic amplitudes over time.

%% Generate a signal with harmonics.
% Build fundamental and harmonic signals.
samplePeriod = 0.001;
time = 0:samplePeriod:5; % [s]
harmonics = [1, 3, 5, 7];
fundamentalFreq = 50; % [Hz]
componentFreqs = fundamentalFreq .* harmonics; % [Hz]
componentFreqs = componentFreqs .* 2 .* pi; % [rad/s]
componentFreqs = repmat(componentFreqs', 1, numel(time));
componentSignals = sin(componentFreqs .* repmat(time, numel(harmonics), 1));

% Apply decay to the harmonics.
decay = [0, 0.1, 0.5, 1];
harmonicAmplitude = zeros(numel(harmonics), numel(time));
for iHarmonic = 1 : numel(harmonics)
    harmonicAmplitude(iHarmonic, :) = exp(-1 * decay(iHarmonic) * time);
    componentSignals(iHarmonic, :) = componentSignals(iHarmonic, :) .* harmonicAmplitude(iHarmonic, :) .* 1;
end

% Generate final signal.
signal = sum(componentSignals, 1);

%% Perform time-frequency analysis.
sampleFrequency = 1 / samplePeriod;

% Compute the spectrogram.
NFFT = 2^12;
winDuration = 0.1; % [s]
winLength = winDuration * sampleFrequency;
winFactor = mean(hamming(winLength));
[S, freq, specTime, p] = spectrogram(signal, winLength, [], NFFT, sampleFrequency);

% Compute FFT amplitude.
amplitude = abs(S) .* 2 ./ winLength ./ winFactor;

% Threshold.
floorLevel = max(max(amplitude)) * 0.1;
amplitude(amplitude < floorLevel) = 0;

% Determine the amplitude of the harmonics of interest.
fundamentalFreq = 50; % [Repeat of above]
harmonicsOfInterest = 1 : 10;
harmonicSearchWidth = 5; % [Hz]
harmonicBracketLower = harmonicsOfInterest .* fundamentalFreq - harmonicSearchWidth / 2;
harmonicBracketUpper = harmonicsOfInterest .* fundamentalFreq + harmonicSearchWidth / 2;

% Find the amplitude maxima for each harmonic at each point in time.
estimatedHarmonicAmplitude = zeros(numel(specTime), numel(harmonicsOfInterest));
for iHarmonic = 1 : numel(harmonicsOfInterest)
    activeBracket = freq >= harmonicBracketLower(iHarmonic) & freq <= harmonicBracketUpper(iHarmonic);
    bracketAmplitude = amplitude(activeBracket, :);
    estimatedHarmonicAmplitude(:, iHarmonic) = max(bracketAmplitude, [], 1)';
end

% Visualise output.
figure;
surf(specTime, freq, amplitude, 'EdgeColor', 'none');
xlabel('Time [s]');
ylabel('Frequency [Hz]');
zlabel('Amplitude');
view(2);

figure;
plot(specTime, estimatedHarmonicAmplitude);
xlabel('Time [s]');
ylabel('Amplitude');
box on;
grid on;
legendBuilder = @(a) ['Harmonic: ', num2str(a, '%1.0f')];
for iHarmonic = 1 : numel(harmonicsOfInterest)
    legendString{iHarmonic} = legendBuilder(harmonicsOfInterest(iHarmonic));
end
legend(legendString);


pairwise error probability and union bound of space time code i have a question that might seem silly.
For a given space time code, say a constellation of 16 unitary matrices. When using the union bound to derive an upper bound on the average error probability (using the equation in the image below)

is it logical to have the Bit Error Rate result higher than 1 for a given SNR? even if what i am deriving is an upper bound ?

Thank you

The union bound can indeed result in a BER higher than 1. As you say, it is just an upper bound, and it becomes tighter as the SNR increases. For low SNR, the union bound can be quite loose.


upsampling for a signal We have a signal $s(t)$.

If we do an upsampling, does the signal duration increase?
What is the point of upsampling if the signal time increases?
Can we do an upsampling if we don't use a shaping filter?

No, if you upsample a signal, its duration (i.e. measured in real-world units like "seconds") does not change.

If your signal $s[n]$ is sampled at a rate $f_s$ Hz, and it is $N$ samples long, then the total duration of the signal is $\frac{N}{f_s}$ seconds.
Say you upsample the signal by a factor $M$, yielding a new signal $s_u[n]$. It will be $MN$ samples long, but it will be sampled at a new rate $Mf_s$. Its duration is therefore $\frac{MN}{Mf_s} = \frac{N}{f_s}$. You can see that the signal duration is the same as $s[n]$.

It's not clear what you mean by a "shaping filter" in your third question.


FFT of random binary data I am trying to make sense of FFTs and binary data.

Say I have a series of random binary data, which is measured with a repetition rate of 400Hz (interval time of 0.0025s). I have a total of 12489 points, which corresponds to a total measurement time of about 31 seconds.

I would like to be able to learn more about what I would expect an FFT of this data to look like.

Some things I would like to understand the significance of:

What should be the average amplitude of the data, post-FFT?
what is the significance of the maximum amplitude of binary data that is not random, but consists of 1,0,1,0,1,0 data (12489 points). How can this help me find my y-scale?

I had a similar question myself, so I created a html demo to test this.

I think you may find it useful in understanding what the FFT is actually doing since you can also read and edit the source code.

https://seanwasere.gitlab.io/FFTBinary

This demo is doing FFT analysis of binary data. 1s and 0s

In my demo I have 1024 bits, and I convert them into an 8 bit wav, and then analyse it using the browsers inbuilt audio api fft analysis function. Firefox browser works best in my case.

The 8 bit wav is

1 channel
samplerate 48000
1024 bits = 0.021333333333333332 seconds in length
the 1s and 0s are actually converted behind the scenes to 256s and 128s

In the user interface there are 4 panels

1st is the FFT result of the binary data drawn as a spectrum graph,
2nd is the FFT result drawn as a waveform. It is essentially triangular spikes or squares, depending on the sequence 0,1,0 or 0,1,1,0.
3rd panel is a view of the actual binary data that was converted into the waveform before FFT analysis.
4th panel is your options, you can

set to have a '1' every N bits by using modulus,
and/or play random 1s and 0s which I call binary white noise.
set the FFT buffer size, the lower numbers such as 128 give lower quality output but is much faster to compute.

The source code is MIT so you can play around with it.


How do I interpret an image of a spectrogram of an audio wind signal? I am running the following code in python to get the spectrogram of a audio signal of wind:
import librosa
import numpy as np
#  Load the audio as a waveform `y`
#    Store the sampling rate as `sr`
file_path = "wind.wav"
y, sr = librosa.load(file_path, sr = None, mono=True, offset = 0.0,  duration=None)
D = librosa.stft(y)
librosa.display.specshow(librosa.amplitude_to_db(D,ref=np.max),y_axis='log', x_axis='time')
plt.title('Power spectrogram')
plt.colorbar(format='%+2.0f dB')
plt.tight_layout()

I get the following image:

My question: How do I interpret this image? What exactly am I looking at and what are the frequencies in this audio file?

The almost steady colored part between 1 s and 35 s, and below about 30 Hz may be the sign of a   low-frequency, continuous noise that seems to be common to mild winds. There seems to be occasional gusts of wind around 12 s to 17 s or 22 s to 25 s. However, the needle aspect of the  spectrogram coefficients may suggest to modify the window, and the little information on frequencies above 256/512 Hz would suggest to look closer at the low frequencies, possibly with a linear frequency scale to better see details.


mathematical expression to detect modulated data within a vector Assume I have such modulated data which is, for example, $x=0.7+0.7i$. That modulated data is encapsulated in a vector as below:

OR

where $c$ is any constant number, let’s say, for example:  $c=0.7$ . My question is how can I detect if the received vector is $y_1$   or $y_2$ ? I mean is there a mathematical expression I can apply for detection?

The same as you'd do for a scalar reception:

Build a stochastic model for your transmit vectors
Build a stochastic model for your noise
from 1. and 2., build a stochastic model for your received vectors
Define a metric on the error statistics that you want to minimize
from 3. and 4., derive a decider

For example:

$x$ are points from an amplitude-phase shift keying modulation $\mathcal X$, with two different amplitudes, and 4 points on the inner "circle", as well as 6 on the outer, and the inner points being chosen twice as often as the outer points; points on the same amplitude circle are equiprobable.
From that follows a statistic for $\vec y_1=(y_{11},y_{12})^T$, call it $f_{y_{11},y_{12}}(\xi_1, \xi_2)$.
In many cases this is relatively easy, since the possible points are discrete, in others you'll have to employ more complicated models e.g through random variable pdf transformations using the Jacobi determinant.
Noise $\vec n = (n_1,n_2)^T$ is uniform in phase over $[0,\pi/2]$, and exponentially distributed in amplitude, strongly stationary, white and independent from $\vec y$. That means it can be described by a single $f_n(\chi)$.
$\vec r = (y_{11}+n_1, y_{12}+n_2)^T$. If you modeled your noise to be independent from signal in 2., then the noise pdf is just the convolution of $f_{\vec y}$ and $f_{n}$.
You define that you want to minimize the symbol error probability for every receive point. Other choices might include minimizing bit error rate locally, or symbol error rate globally under constraints (e.g. complexity).
You hence decide to build an MAP decider that gives you an estimate $\hat x$ that gives you $\hat x_{\vec r} = \operatorname*{\arg\max}\limits_x f(\vec r|x)$. So, you write down the 2D density, and find a maximum.


How can I convert a step response that is described by the sum of two exponents to an IIR filter analytically? If a system has the following step response:

$$\begin{align}
y(t) &= 1 - \Big(A e^{−t/T_1} + (1-A) e^{−t/T_2}\Big) \\
 &= 1 - A e^{−t/T_1} - (1-A) e^{−t/T_2} \\
\end{align}$$

and I know the values of $A$, $T_1$ and $T_2$, can I transform this to a discrete-time IIR filter analytically?  If there is no analytical solution, is it possible to say how many poles/zeros the IIR filter should have?

I've done it numerically by using the Prony's method on the impulse response, but I'm hoping there is a more direct way.

You can use the step-invariant transformation, which makes sure that the step response of the discrete-time system matches the step response of the continuous-time system at the sampling instants. I suppose this is what you mean by "exact solution".

Sampling the given continuous-time step response at $t=nT$ gives the desired step response of the discrete-time system:

$$g[n]=y(nT)=\left[1-Ae^{-nT/T_1}-(1-A)e^{-nT/T_2}\right]u[n]\tag{1}$$

where $u[n]$ is the unit step.
The $\mathcal{Z}$-transform of $(1)$ is

$$G(z)=\frac{1}{1-z^{-1}}-\frac{A}{1-e^{-T/T_1}z^{-1}}-\frac{1-A}{1-e^{-T/T_2}z^{-1}}\tag{2}$$

The frequency response of the desired discrete-time system is then given by

$$H(z)=(1-z^{-1})G(z)\tag{3}$$

After performing some basic algebra and using $a=e^{-T/T_1}$ and $b=e^{-T/T_2}$ we obtain

$$H(z)=\frac{[1-b+A(b-a)]z^{-1}+[a(b-1)+A(a-b)]z^{-2}}{1-(a+b)z^{-1}+abz^{-2}}\tag{4}$$

As a result, the IIR filter with transfer function $(4)$ has a step response that equals the given continuous-time step response $y(t)$ at the sample instants $t=nT$:


Background Noise Reduction in Audio signal I using high pass filter  for reduction or removal of noise of audio signal but still some noise is there in signal.i have MIC,Head Set,Bluetooth and phone are connected to controller,while talking on MIC some background is there.How to reduce the noise and is there any algorithm for it.let me know please.

Im implementing in c code structure.

Any Algorithms ?

Depends a bit on the type of noise. If it's hum or buzz, you need to clean up your recording electronics and if push comes to shove you may be able to notch it out.

Real acoustic background noise is difficult to remove. You may want to read up on "blind source separation". Here are a few examples http://cnl.salk.edu/~tewon/Blind/blind_audio.html


Power in a frequency band: Spectral estimation versus time domain I wish to find total power in a given frequency band within different time windows of a signal.
I see two methods of doing so, first in the frequency domain. One could estimate the power spectral density within the windows (i.e. tapering and then apply FFT), and integrate over the frequency band of interest.
In the time domain, one could apply a bandpass filter to the whole signal, and then estimate power as the mean square of the signal within the windows.
What would be the differences between the results of the two methods, could one of them be considered superior in this case?

Personally I'd go with the bandpass filter, unless there is a strong reason not to.

Much lower implementation complexity, CPU and memory usage.
You can get a continuous output signal at the same sample rate as the input signal. The FFT method will get you only one output per frame, so it's down sampled to the hop size.
No framing artifacts based on how your signal lines up with the frame grid: You can frame it, but you don't have to.
You can design fairly flexible "ballistics" into the detector: fast attack, slow decay, peak hold, lossy peak. I.e. it's easy to customize how fast/slow it reacts to different types of energy changes.

On the downside the choices for frequency weighting functions are fairly limited. If you need fancy frequency domain shapes the FFT may be better.


What is the role of complex exponential? What is the role of complex exponential $ e^{jθ} $ in Fourier Transform? Is it different in the continuous and in discrete time domain?

Euler's relationship says that $e^{j\Theta}$ is equal to $cos(\Theta) + j*sin(\Theta)$.  The Fourier Transform can then be seen as correlating the signal with sinusoids at various frequencies.  The continuous Fourier Transform correlates with an infinite number of sinusoids, while the discrete transform uses $N$ sinusoids, where $N$ is the length of the transform.


Kalman filter for tracking obstacles I'm developing a Computer Vision project in Matlab as an aid to visually impaired people. The setup would be a stereo pair of cameras which the blind person would carry. Using this stereo information, I generate a disparity image using Semi-Global Block Matching (SGBM). After that, I process a "virtual disparity", which basically is a homogeneous transformation, to pass the image plane to the ground, resulting in an image like the following:

Now, After removing the ground plane, morphologically improve the image, and threshold it to keep only near obstacles, I get a foreground mask with the approximate shape of the objects.

I was thinking on tracking each one of the obstacles to improve the detection. Considering that the camera is not fixed since the blind person carrying it is moving, I wonder if the relative movement of all the obstacles in the image would be easily followed by a Kalman filter, or it would be better a feature-based tracker like Kanade–Lucas–Tomasi feature tracker (KLT).

Thank you in advance

You might consider a particle filter.  Here's a link to a paper I wrote about tracking objects in video using a particle filter.  The great thing about these is that objects can be tracked through temporary occlusions.  The trick with using a Kalman filter here is dealing with the nonlinearity introduced by the edges of the video field, and casting your measurement model as a linear combination of your states (presumably x,y location).  The Unscented Kalman filter is another route you could go, it doesn't require linear models.


Which of DFT, DCT and DWT transforms is more robust to noise and geometric transformations in image processing I am currently studying about image watermarking and I have been testing these 3 domains to embed the watermark at.

My question is, how do these domains hold up to noise addition or geometric transformations, regardless of watermark embeding.

To be more precise, I want to find the relationship between an original image and a transformed one (with the transformations being the ones I mentioned in the title, where geometric transformations include rigid and non-rigid ones) in these domains, and which one holds best at minor and/or major changes.

Also, I would like to learn about the generic pros and cons of these domains, compared to one another, aswell as why it is said to be better to use one instead of the others in watermarking.

Sources would be apreciated. Thank you!

From the perspective of numerical analysis, the most robust transformations are unitary orthogonal matrices because they have a condition number of one.  Indeed the condition number is the metric that defines robustness.  You can find this in Stuart’s Introduction to Matrix Computations or the more advanced book by Golub and Van Loan.

A related way of looking at perturbation sensitivity are the Holder p norms for Matrices also covered in Stewart


Zero peaks in fft data even after subtracting mean I have a time series sensor data. It is for a period of 24 hours sampled every minute (so in total 1440 data points per day). I did a fft on this to see what are the dominant frequencies. But what I got is a very noisy fft and a strong peak at zero.
1. I have already subtracted the mean to remove for the DC component at bin 0. But I still get a strong peak at zero. I'm not able to figure what could be the other reason or what should I try next to remove this.
2. The graph is very different from I have usually seen online as I was learning about fft. In the sense, I'm not able to see dominant peaks like how it is usually seen. Is my fft wrong? Attaching code that i tried and images:

import numpy as np
from matplotlib import pyplot as plt
from scipy.fftpack import fft,fftfreq
x=np.linspace(29,32,1440).tolist()
x=np.array(x)
x=x-x.mean()
N = 1440
# sample spacing
T = 1.0 / 60
yf = fft(x)
yf_abs = abs(yf).tolist()
plt.plot(abs(yf))
plt.show()
freqs = fftfreq(len(x), 60)
plt.plot(freqs,yf_abs)
plt.show()

Since I'm new to this, I'm not able to figure out where I'm wrong or interpret the results. Any help will be appreciated. Thanks! :)

You must be assuming you have a dominant peak at zero from only looking at the same graph you shared, but if you really did remove the mean, then value at bin 0 will be 0 (as bin 0 is directly proportional to the mean). Inspect the data carefully as their does not appear to be anything wrong with the code with respect to the FFT.

What is likely occurring is the bins close to bin 0 are non-zero, and these represent the lowest frequency fluctuations in your data. Each of these is zero-mean, but represent the slowest changing frequency components. In this case you are getting the exact results expected, representing what is in your actual data.


Strange behavior from frequency to angular frequency in FFT I have written a code to compute the Fast Fourier Transform of a simple complex exponential with frequency $f=50.0$, using scipy.fft. The code is written below:
import numpy as np
from scipy.fft import fft, fftfreq, fftshift
import matplotlib.pyplot as plt

T = 10
dt = 0.001
f = 50.0

x = np.arange(0.0, T, dt)
y = np.exp(1.j * 2 * np.pi * f * x)

xf = fftfreq(len(x), dt)
xf = fftshift(xf)

yf = fft(y)
yf = fftshift(yf)

plt.plot(xf, 1.0/len(x) * yf, label='yf')
plt.plot(xf, 1.0/len(x) * np.abs(yf), label='abs(yf)')
plt.xlabel('f')
plt.legend()
plt.show()

The output of the calculation is a frequency spectrum with a single peak, for both yf and its absolute value (which are perfectly overlapped) centered at frequency $f=50.0$:

Now instead of using the frequency, I assume to work with an angular frequency $\omega=50.0$:
import numpy as np
from scipy.fft import fft, fftfreq, fftshift
import matplotlib.pyplot as plt

T = 10
dt = 0.001
omega = 50.0

x = np.arange(0.0, T, dt)
y = np.exp(1.j * omega * x)

xf = fftfreq(len(x), dt)
xf = fftshift(xf)

yf = fft(y)
yf = fftshift(yf)

plt.plot(2 * np.pi * xf, 1.0/len(x) * yf)
plt.plot(2 * np.pi * xf, 1.0/len(x) * np.abs(yf))
plt.show()

In the plot I have multiplied xf by $2 \pi$ to convert from frequency $f$ to angular frequency $\omega$.
Surely this corresponds to a frequency $f = \frac{\omega}{2 \pi}$ different from the first example but I would have expected similar behavior. In this case yt have a kind of dispersive shape, while its absolute value is very broad, it is not normalized and it is not perfectly peaked at $\omega=50.0$ but it is slightly shifted to $50.26$.

Could someone give me some explanation about this unexpected behavior?
If I have to work with raw data, coming from a complex calculation, expressed in terms of angular frequency $\omega$, how can I obtain trusty results from the second programs?

This is probably the single most asked question on this forum. A similar one was asked  and answered just a few hours ago: Spectral leakage from mathematical point of view
Your expectations are wrong: you only get a spectral dirac impulse if the frequency of the complex exponential on he FFT grid, i.e. an integer multiple of the sample rate divided by the FFT length. In your case you have $N_{FFT} = 10000$ and $F_s = 1000$, so your FFT bin spacing is 0.1 Hz. 50 Hz is an integer multiple of 0.1Hz but $50/(2\pi)$ isn't and so the energy will spread out over the adjacent bins.
A different way to look at it: the DFT assumes that the time signal is periodic with the FFT length. Since your exponential does NOT have an integer number of periods inside one FFT frame, the periodic repetition creates a discontinuity at the frame boundary


Minimum Energy Delay Property I have 4 systems which are represented in z domain like this:

Where I know the poles and zeros.

Then I compute Partial Energy Sequence of impulse response of all these system using this formula:

And I get this result:

And I want to observe minimum energy delay property of the minimum phase systems can be specified as:

If I am right, the formula above tells us that minimum phase system has a bigger or equal at least partial energy but, unfortunately, you can see in partial energy sequence plot, the minimum phase system has the smallest partial energy!

Am I interpreting the minimum energy delay formula wrong or what am I missing here?

Given only the zeros and poles, the impulse response is determined except for a constant $K_i$, which you don't show.

From your calculations, it is clear that the energy of the impulse response for all filters is different.

If you normalize all filters so that all impulse responses have the same energy (i.e. $\sum_n |h_i[n]|^2$ is the same for the 4 filters), then the 4 curves for partial energies will end up at the same value, and it will be clear that $H_1$ has its energy distributed earlier in time than the others.


HEVC Deblocking filter algorithm I'm doing research on block noise caused by video compression (especially in newer codecs) and I've read the IEEE article about the HEVC deblocking filter (here).

However, when reading the algorithm I'm confused because the formula seems wrong. The first formula to determine if the deblocking filter will affect a block is:

$ |p_2,_0 - 2p_1,_0+p_0,_0|+|p_2,_3 - 2p_1,_3+p_0,_3|+|q_2,_0 - 2q_1,_0+q_0,_0|+|q_2,_3 - 2q_1,_3+q_0,_3|>\beta \ (1)$

It is also written (and shown in a figure) that $\beta $ increases with the quantification parameter QP.

Moreover, later it is written that to choose between normal and strong deblocking, the following equation is used (with $i$ going from $0$ to $3$):

$ |p_2,_i - 2p_1,_i+p_0,_i|+|q_2,_i - 2q_1,_i+q_0,_i|<\beta/8 \ (2) $

Since IEEE is a respectable source and the paper is a couple years old, it would surprise me that there is a typo. However, I do not understand how this formula can work. It does seem to me that the second inequation would always be wrong at least for some $i$ if the first one is true.

Am I misunderstanding how the filter works or is it a typo in the paper?

The ffmpeg implementation of the HEVC decoding follows this formula instead:

$ |p_2,_0 - 2p_1,_0+p_0,_0|+|p_2,_3 - 2p_1,_3+p_0,_3|+|q_2,_0 - 2q_1,_0+q_0,_0|+|q_2,_3 - 2q_1,_3+q_0,_3|<\beta \ (1)$

This appears in the hevc_deblock.asm file.

;compare
    pcmpgtw         m15, m13, m14
    movmskps        r13, m15 ;filtering mask 0d0 + 0d3 < beta0 (bit 2 or 3) , 1d0 + 1d3 < beta1 (bit 0 or 1)
    test            r13, r13
    je              .bypassluma

So I conclude that there is a mistake in the paper and the correct formula is that one (since it also makes more sense).

edit: As the comment suggested, I checked in the reference implementation for the algorithm and the interesting part is here in the file TComLoopFilter.cpp:

for (UInt iBlkIdx = 0; iBlkIdx<uiBlocksInPart; iBlkIdx ++)
      {
        Int dp0 = xCalcDP( piTmpSrc+iSrcStep*(iIdx*uiPelsInPart+iBlkIdx*4+0), iOffset);
        Int dq0 = xCalcDQ( piTmpSrc+iSrcStep*(iIdx*uiPelsInPart+iBlkIdx*4+0), iOffset);
        Int dp3 = xCalcDP( piTmpSrc+iSrcStep*(iIdx*uiPelsInPart+iBlkIdx*4+3), iOffset);
        Int dq3 = xCalcDQ( piTmpSrc+iSrcStep*(iIdx*uiPelsInPart+iBlkIdx*4+3), iOffset);
        Int d0 = dp0 + dq0;
        Int d3 = dp3 + dq3;

        Int dp = dp0 + dp3;
        Int dq = dq0 + dq3;
        Int d =  d0 + d3;

//...

        if (d < iBeta) //<-not a >Beta as the paper said
        {
          Bool bFilterP = (dp < iSideThreshold);
          Bool bFilterQ = (dq < iSideThreshold);

          Bool sw =  xUseStrongFiltering( iOffset, 2*d0, iBeta, iTc, piTmpSrc+iSrcStep*(iIdx*uiPelsInPart+iBlkIdx*4+0))
          && xUseStrongFiltering( iOffset, 2*d3, iBeta, iTc, piTmpSrc+iSrcStep*(iIdx*uiPelsInPart+iBlkIdx*4+3));

          for ( Int i = 0; i < DEBLOCK_SMALLEST_BLOCK/2; i++)
          {
            xPelFilterLuma( piTmpSrc+iSrcStep*(iIdx*uiPelsInPart+iBlkIdx*4+i), iOffset, iTc, sw, bPartPNoFilter, bPartQNoFilter, iThrCut, bFilterP, bFilterQ, bitDepthLuma);
          }
        }


Complex FFT reversing I've succeeded to compute complex FFT of a data array of interleaved values using the arm_cfft_f32() API of the CMSIS DSP library for Cortex-M4, and got in the same array the FFT results as frequency domain values with the same interleaving.

I chose a result value of the Frequency Domain (let say a point) based on a specific criteria, and want to go back to the initial index (in the Time Domain) before calculating the FFT and passing through arm_cfft_radix8by4_f32() API.

Is there anyway to go in the opposite sens, and identify a point based on its FFT result value ?

Thanks and regards,
tinymanta.

I think you're misunderstanding the relationship between the time and frequency domains. Each point in the output of a DFT is a function of all of the input time-domain values. There is no one-to-one mapping between input points and output points.


Does the DFT calculate spectral components up to half the sampling frequency, $f_s/2$? This question is prompted by a statement made in this response (reproduced below):

  The DFT calculates spectral components up to $f_s/2$, no matter what the input signal is.

A book I'm reading explains that the DFT of a signal $x(n)$ is the result of the inner product of $x(n)$ with each of the sampled complex sinusoidal basis vectors $s_k(n)$. Here is the author's definition of the DFT:

The value $k$ is called the bin number and $\omega_k$ are harmonics of the sampling frequency $f_s$ (expressed in angular unit of radians, alternatively represented as $f_k = \frac{k}{N}f_s$).

Since $k$ is taken from $0$ to $N-1$, $X(\omega_k)$ is defined up to and including $f_{N-1} = \frac{N-1}{N}f_s$. Why would Deve (answer author) say $f_s/2?$

Depends a bit on the indexing convention, but in typically you would interpret the frequency interval as $[-f_s/2, f_s/2]$ and not as $[0,f_s]$. The DFT is periodic in N so you you have

$$X(N-1) = X(-1) $$

Keep in mind that that the sampling theorem requires the signal to be band limited to $f_s/2$ so assuming that you have actual independent information about frequencies higher than $f_s/2$ is misleading.

For real signals, it's conjugate symmetric anyway, i.e. $f_{-k} = f^*_k$ so there is only independent information on half of the spectrum anyway.


GNU Radio GMSK modulator I'm looking at the GNU Radio GMSK modulator block. The Gaussian taps are somehow convolved with a rectangular window with 4 * samples-per-symbol taps owing to the fact for GMSK, the filter spans four symbols. I was wondering if anyone could elaborate why do we convolve the Gaussian taps with the square wave?

self.gaussian_taps = filter.firdes.gaussian(
        1,             # gain
        samples_per_symbol,    # symbol_rate
        bt,            # bandwidth * symbol time
        ntaps                  # number of taps
        )
self.sqwave = (1,) * samples_per_symbol       # rectangular window
self.taps = numpy.convolve(numpy.array(self.gaussian_taps),numpy.array(self.sqwave))
self.gaussian_filter = filter.interp_fir_filter_fff(samples_per_symbol, self.taps)

Let's say you want to filter a signal $x[n]$ through a Gaussian filter with impulse response $h_g[n]$ and a Moving Average (the "sqwave") filter with impulse response $h_s[n]$.  Then the resulting operation is

$$y[n] = h_s[n] * (h_g[n] * x[n]) = (h_s[n] * h_g[n]) * x[n]$$

So the a priori convolution of the two filter tap sets into one filter taps set, just reduces the two filtering operations into one filtering operation.

The Moving Average filter is just a low pass filter with a poor transition band and poor sidelobe supression.  (The one used here isn't even normalized to unity gain.)

I can't definitely answer "why?", except to say that it looks like someone wanted to add a cheap low pass filter as a sort of channel filter.  That moving average filter has been there since 2007, when GNURadio's git history starts.  The old CVS repositories aren't on-line anymore, so the original rationale is probably lost to history.


What is the ROC for this discrete signal: $$ x(k)=4[u(k-2)-u(k)*δ(k-3)]$$

I found that the $\mathcal{Z}$ transform of the signal is $X(z)=4/(z^2)$.

What would the ROC be?

HINT:

Remember that the ROC is the region in the $z$-plane for which the series

$$X(z)=\sum_{n=-\infty}^{\infty}x[n]z^{-n}\tag{1}$$

converges. Since you've found $X(z)$ you also know that $x[n]=4\delta[n-2]$, i.e. there's only one value of $n$ for which $x[n]$ is not equal to zero. What does that mean for the convergence of $(1)$?


Obtaining the variance of the image noise I am implementing a digital image correlation algorithm, and I found that for some calculations the variance of the image noise would be very interesting.

How can I determine that parameter? The only thing I have are a set of 1000 images of a metal deformation process and no prior information about the noise in the image.

I found the answer.

The correct approach is the next one:

In a controlled environment (appropiate ligth, no vibrations,...) take several images of the same object (e.g. a white sheet of paper).
Then calculate the variance of each pixel on time (in different images).

Finally mean all those variances and there it is! the variance of the image noise.


What exactly does compression say about correlation of data? I've been using the following formula on various empirical data $d$, to obtain a correlation factor $c_f$:-

$$ c_f = { |C(d_s)| \over |C(d)|} $$

where $C$ is a compression function like bz2 or zip, and $d_s$ is the data having been randomly shuffled. The purpose of the random shuffle it to destroy any byte wise auto correlation that might exist in $d$.

A typical value for $c_f$ might be 2.7.  What might this indicate viz a viz any possible correlation within the data?  Is it at all possible to convert 2.7 to a lag figure of $n$ bits, even loosely?

Example figures are:-

$|d|$ = 921,600 bytes

$|C(d)|$ = 201,000 bytes

$|C(d_s)|$ = 543,000 bytes

What might this indicate viz a viz any possible correlation within the data?  Is it at all possible to convert 2.7 to a lag figure of n bits, even loosely?

The compression ratio is directly proportional to the entropy ($H$) of the source. For example, Shannon's Entropy tells you how many bits on average would be required to encode a source, given its statistical distribution. The expression for Shannon Entropy is:

$$ H(S) = -\sum_{i=0}^{|S|-1} P(s_i) \log_2 P(s_i)$$

Where $|\cdot |$ indicates the size of an ordered set, $S$ is a memoryless source having $|S|$ symbols, each one appearing with probability $P(s_i)$ in a stream.

If the source was totally predictable ($1111111111111111111111111111 \ldots$) then $S_p=[1]$, $P(S_p) = [1]$ and to describe it you would need $H(S_p)=0$ bits.

If the source had to alternate between two states but in a totally unpredictable way then $S_u=[0,1], P(S_u) = [0.5,0.5]$ and to describe it you would need $H(S_u) = -(0.5 \cdot \log_2 0.5) \cdot 2 = 1$ bit.

BUT, if the same source started showing a preference for one of the states, it suddenly became just a bit more predictable. So, let's say $S_b=[0,1], P(S_b)=[0.7,0.3]$ and $H(S) = -(0.7 \cdot \log_2 0.7 + 0.3 \cdot \log_2 0.3 \approx 0.88129...$.

Remember that $H$ is expressed in bits here. If I need $H(S)$ bits per symbol in $S$ and produce a message $M$ that is $N$ symbols long, then I need $|M| = H(S) \cdot N$ bits to describe the message.

The unpredictable source $S_u$ needs 1 bit per symbol, so, we need as many bits as symbols in $M$ and $|M_u| = N_u \cdot H(S_u)$. But the "biased" source $S_b$ needs $\approx 0.88129$ per symbol, so $|M_b| = N_b \cdot H(S_b)$ and theoretically at least, we could fit it in less space.

The question now is, how do you do that? How do you re-encode a sequence to occupy less space? Because that is very, very important (for the interpretation of metrics such as this one).

Without wanting to go into too much detail (now), let's say that we do have a way of encoding a message in such a way that the compression ratio approximates its entropy as the length of the message tends to infinity. This is what you call $C$ (and it is encouraging that you have selected zip, bz2, etc).

Then, to form a ratio of the lengths of the shuffled and unshuffled compressed sequences is equivalent to forming the ratio of $\frac{|M_{u \leftarrow b}|}{|M_b|}$.

In other words, you have a file that is coherent because it bears content. Because of this reason, the bytes appear in some predictable way, following some kind of pattern, kind of like $S_b$ but this time around your $S_b$ contains 256 symbols and your $P(S_b)$ contains 256 probabilities (this is not exactly accurate here because a file is not memoryless but I will come back to this later).

When you shuffle the file, you basically construct an equivalent white noise sequence from a source $S_u$ where $P(S_u) = \frac{1}{256} \forall s \in S_u$ (Where $\forall$ means "for all"). This is what I mean up there with that $M_u \leftarrow M_b$.

Earlier, $S_u,S_b$ were unrelated, but here $S_u$ is the result of $S_b$ after you shuffled it.

So, having said all this, what you basically have expressed there is a ratio that is equivalent to the inverse of the compression ratio of $C$ on your data $d$.

Provided that the shuffling is adequate $|C(ds)|$ is equivalent to the length of $ds$ because it is impossible to compress it. And, the length of $|C(d)|$ is approximated by the entropy of $d$ times its length.

What does this tell you? It tells you how much bigger is a white noise sequence using the same symbols in $d$ compared to the coherent version of $d$.

What does this say about the auto-correlation of the data? Under very strict conditions it can tell you the average size of the longest redundant sequence of bits. If the data compression ratio is 0.5, on average, you are substituting 2 bits for 1.

And this, is where things can go really bad, in terms of choice of $C$ and interpretation.

The compressor $C$, doesn't care what is in $d$. The only thing it sees is a stream of bytes. If you have a file of doubles in a common format, the compressor will compress each and every one of these 8 bytes separately, not as one number. What your ratio will return, is information on the efficiency of compressing the stream, not your "data".

  Is it at all possible to convert 2.7 to a lag figure of n bits, even loosely?

This was answered above, but this discussion so far brings us to the wonderful subject of LZ Complexity.

The LZW algorithm is a dictionary based algorithm. It compresses a sequence by keeping a running dictionary of fragments it has encountered before.

Consider 0101001010010. We start at the first 0, have we seen it before? No, put it in the dictionary, output 0, move to the next. Have we seen 1 before?, No, put it in the dictionary, output 1, move to the next. Have we seen 0 before? Yes, carry on to the next. Have we seen 01 before? No, put it in the dictionary, output the index of the known sequence (the 0) and the new bit (1) and move to the next....and so on until the end of the file.

You see what's happening here, the algorithm compresses when it encounters a sequence it has seen before and it substitutes it for its index in the dictionary. In these early stages it does not seem to compress much but if the source is predictable then the algorithm will build a set of increasingly long repeating sequences very quickly and it will be able to start substituting repetitions. Have we seen 001100110011001100110011 before? Yes, it is entry 242. So, you substitute a long string of bytes for an integer.

The reason I am talking about this is because if you wanted to make such strong statements about the lag of the autocorrelation, it is not enough to examine the compression ratio. Ideally, you would have to examine the patterns within the dictionaries.

If you are lucky, your data is inherently digital (categorical) and you might be able to do this. But if your data is real numbers, then you have to quantise them and the way you do this makes it very difficult to apply these methods to real valued data.

Hope this helps.


Time-frequency analysis of non-sinusoidal periodic signals Given the history of the sum of a time-varying mixture of periodic signals, say square waves, how would you efficiently estimate the number and frequencies of components active at a particular time? The amplitudes and frequencies of the components are arbitrary but fixed real numbers; if a component is active at a certain moment, it will retain its amplitude and frequency if it is active again.

Let us start with the unsupervised methods...

A first approach would be to compute a spectrogram and factorize it with NMF (Non-negative Matrix Factorization). If you are unfamiliar with this technique, it decomposes a spectrogram into a sum of $k$ constant-spectrum sources, each of these having a time-varying amplitude envelope applied to them. This model perfectly suits your problem and there is a good chance that the columns of the decomposition will be spectra for each  waveform/frequency pair (from which you can estimate a fundamental frequency through spectral sum or by turning it back into an autocorrelation function), and the rows will be the activation signals. There are many implementation out there and it's cheap to just throw your problem at it, there's a good chance it'll just work.

Note that the underlying model of NMF is less restrictive than your signal model...

First, because your activation envelopes are either 0 or 1. I haven't dealt with such constraints in the past, but there's probably a way of defining a penalty measure on the terms of the activation matrix $W$ different from 0 and 1, adding that to the optimized criterion and derive a new set of multiplicative update equations. One constraint not mentioned in your comment is that you probably also assume that the sources are not "blinking" in and out rapidly, and stay active/inactive, for a rather long number of consecutive frames. Constraints can be added to NMF to penalize discontinuous activation envelopes. See Virtanen's paper on continuity constraints for that.

The second specificity of your problem is that you might want to make the assumption that all sources have the same waveform, that is to say, that all sources have the same log-frequency spectrum, modulo a translation. To address this specifity, a recommended technique is to compute a constant-Q spectrogram on your input data (which has the right shift-invariance), and perform a shifted-NMF, with one target source. This will recover the spectrum common to all your sources; and the spectral shifts and temporal activation matrix will provide the equivalent of a spectrogram tailored to your source spectrum. See the synthetic signal in the results section of this paper - this is exactly the signal model you describe. This might require some work, but you could probably tweak this method as well to your more rigid signal model in which the activations are either 0 or 1.

Now on to the supervised methods.

If the set of candidate waveforms is known in advance and small and/or the number of frequencies at which they can appear is small (say you want to transcribe music played on a 80s toy keyboard), you could afford matching pursuit. This will correlate your signal with impulse responses consisting of a windowed signal for all searched waveforms/frequencies pair. Computational costs might rapidly make it a bad choice.

Another supervised method is to use the same kind of multiplicative gradient updates as for NMF, but keeping the matrix $W$ containing the spectral templates "locked". That is to say, you build a matrix W with all candidate spectra you want to probe, and perform your multiplicative update to decompose your spectrogram into a product of an activation matrix H and W. Again, update rules tweakable to include the binary activation constraints. See Bertin's work on piano transcription which uses this technique to factorize spectrograms into a sum of synthetic piano note spectrograms (a good start here). This has also been used for drum transcription by Paulus - using drum samples as a basis.


Segmenting page with dense input and disconnected separators I would like to segment pages of a historic book, which is printed densely and shows disconnected black column separators, see below.

What I would like to arrive at is the following segmentation.

I have been playing with ocropy so far, which makes use of scipy/ndimage. But the tool used is not important, I'm interested in methods (morphology, ...) that can solve this.

Any ideas?

1 Top-down approach

Your problem looks like you'd

want to segment horizontally,
shrink your segments from the left and right image border,
check the number of text columns in each horizontal segment, and then
within these segments simply look for the columns.

1.1 – segment horizontally

Laziest approach I could come up with:  simply find completely white consecutive runs of lines.

To avoid detecting the whitespace between normal lines of text, have a "safe" boundary for large inter-paragraph space. Below that, check whether the non-white lines below and above roughly have the same length.

1.2 – shrink horizontally

Simply find the outmost left and right columns that contain no black. Crop these and everything beyond them out.

1.3 & 1.4 – check number of columns & vertically segment

I'd go for: find columns in the segment of pure white, approach "safety margins" just like for the horizontal segmentation. If you're left with a 1-3px wide vertical segment, check whether it roughly looks like a vertical line.

Another approach would be to directly search for vertical black lines (i.e. go through columns, count black pixels. If column more than let's say 40% black –> vertical line), or if you have to derotate the images first, anyway, start your whole process prior to 1 with a Hough transform and look for roughly vertical lines close to the center of the image – these would show up as peaks in the Houghes-typical position/angle plane.

2 Bottom-up approach based on content

Another thing I'd try is to logically:

Detect letters' and vertical lines' outlines
Detect the text lines' outlines (by grouping letter outlines and using vertical lines as boundaries)
Detect the paragraphs' outlines (by grouping lines)

2.1 – detect letters

Pretty much: find black pixel, find all directly connected black pixels, find the min and max x and y of these groups, draw red rectangle enclosing them. Repeat for next black pixel.  (with "draw red" I pretty much mean "mark this as already processed")

2.1a – detect vertical lines

The same as 2.1, but simply filter the detected rectangles by aspect ratio. Mark as green, or however.

2.2 – detect lines' outlines

Starting from the left, take a symbol outline, find center, move right till you either

hit a vertical (green) line or
exceed the maximum sensible word distance or
hit the image border or
hit the next symbol.

In cases 1.–3., mark all found symbols as content of one line. In case 4., add found symbol to line "stack", and repeat.

2.3 Find paragraph outlines

From top left, take a line outline. Go down.

If the next outline you hit has roughly the same start and width, add to lines stack. If not, or you hit a boundary, declare all lines on stack a paragraph, mark processed, take next unprocessed line outline from top, and repeat.

3 Bottom-up approach, whitespace-based

Idea: on black and white input

Pick a white pixel
Find largest undisturbed white (or red) rectangle enclosing said pixel. If area smaller than threshold, mark green. If larger than threshold, mark red
Go back to 1.

The result should pretty much be a red mask, that, when inverted, gives you the segmentation you need.


Does receiver always know about the noise level in real world application? In many research articles, the effect of measurement noise on estimation performance is often reported. But it is not clear to me what is the proper way of using the signal to noise ratio (SNR). For example, the variance of measurement noise which is assumed to be Additive White Gaussian Noise( AWGN) is varied to obtain different SNRs. Then for each SNR the mean square error (MSE) between the parameter estimates and the known parameters are calculated. When SNR is high, it means that the amount of noise is less.

The following code in Matlab is based on my understanding that SNR is defined after the input to a system is passed and we get the output from the system, then we add noise for a particular SNR. Also, please correct me if I have put any incorrect information. Thank you.

In the code, I first generate some data $x$ consisting of $N=100$ data points. Then I have randomly generated 3 coefficients representing the impulse response of an FIR system. The data $x$ is passed through the FIR system to obtain $y$. I then add AWGN of SNR = 40dB.

Question1: My question is which step of the estimation stage is the SNR defined?  Can somebody please confirm if this is the correct approach or not?

Question2: If SNR = 40dB, how does one know the variance of the noise at the receiver end? In practice (in industry application) does the receiver end always know about the level of SNR and the variance of the noise?

N=100;
x = randn(1,N); %generating random data
h = rand(1,3); %some unknown parameters representing impulse response
y=filter(h,1,x); %FIR filter
z = awgn(y,40,'measured'); % adding AWGN measurement noise at SNR = 40dB
%Run some estimation method to estimate h_hat

It depends somewhat on the application. For an active sonar system the receiver has to deal with the near backscatter ( reverberation) of the active pulse ( unless it is continuous wave) so there typically is some blanking and fast/slow automatic gain control using feedback.  In a passive sonar receiver, the same philosophy using feedback based AGC is typical. Similar concepts apply to some radar systems.  Down stream, there are a number of ideas used that include leaky integration to track slow variations in background, and if your detection architecture is a large number of cells like fft bins, there are schemes where means or medians of  neighbor cells adjust the gain of a particular cell.

In these kind of systems, the Neyman Pearson  criteria is used to set thresholds, so you really don’t need to know SNR, you need to have your AGC and other schemes give you a known noise variance. While an optimal likelihood ratio would include SNR as a known term, one can go with a locally optimal assumption, like a weak SNR or a generalized likelihood ratio. Neyman Pearson also circumvents having to know signal prior probabilities of occurrence.

Communication receivers also use fast/slow feedback AGC. Systems have a finite dynamic range so you have to put your signal at a reasonable operating point. This isn’t the same as knowing your SNR but an AGC lets you do things like pilot detection. If necessary, SNR becomes an online estimation and tracking problem. A clean signal is often equivalent, in a practical sense, to knowing your noise level.

There is a phenomena known as the “threshold effect” that can be explained as there being a relatively narrow range of SNR increase where performance goes from very poor to very good.  Knowing your exact SNR below this threshold range has no utility, performance is bad regardless. Knowing your exact SNR above this range is probably not consequential either.

The answer to your question is yes, a practical receiver does need to know the noise level, but typically in the context of placing the signal at an operating point compatible with its dynamic range.  The locus of this function is AGC.  There may be additional downstream SNR processing as well. It isn’t necessarily located in a single place.  CFAR is one kind of downstream processing for radar/sonar.


Re sizing and thresholding small font scanned Arabic words I am working on an arabic OCR for printed scanned documents ... Some of the scanned documents may written with font size of height 8 which is quite small ... I want to resize the height to 60 pixels but some artifacts may occur due to nature of arabic characters .. some characters may overlap .I had used local thresholding methods after resizing but the results are still not acceptable ... any ideas ?

This is an example image:

This is the same example after resizing and applying local adaptive thresholding using 50 as a window size:

As you can see there are some discontinuities in some characters like this one :

Is there is any method to resize the image while preserving textual shape ?

My approach to fix characters discontinuities :

Threshold original image before resizing using local adaptive threshold using window size of 16 (this will solve characters discontinuities but the holes in the characters are filled) call it smallbw.
Resize smallbw using imresize(smallbw, [nh nw], 'nearest') and fill the holes in characters using imfill
Resize the original image to height 60 pixels using imresize(originalIm, [nh nw], 'nearest') call it largebw
Fill holes in largebw using imfill and call it bwfill
Extract holes from largebw by bwholes = bwfill - largebw
Finally, subtract bwholes from smallbw to get this

you can see here that the discontinuity found in the character @Image 3 has been solved ... but there is another problem raised here some characters may overlap as shown here

This is the best results I could achieve so far ... are there any other ideas that may solve these issues ??
and if you think this issue has no solution how would I resolve it rather than using resizing? what about using 12 font size text instead of 8?

Useful links :
Local adaptive threshold method used

Operating system : windows 7

Programming environment : Matlab 2013a - Image Processing toolbox

I think you could specify a higher DPI for higher 'resolution' (which effectively changes the pixel number of your image) in loading the data.

Resizing the data will usually cause small changes to your data that are usually irreversible (unless you know exactly the correct parameters for resizing).


Resampling an audio signal I am trying to resample an audio signal in my application. It will be between 48kHz and 96kHz audio (both upsampling and downsampling).

Looking at the Wikipedia page for resampling I see there are two approaches:

Decimation and upsampling
Polynomial interpolation e.g. Lagrange

My questions are:

The upsampling Wikipedia page states

  Smooth out the discontinuities with a lowpass filter, which replaces the zeros.

How does inserting L-1 zeroes between the original sample and then passing it through a low pass filter magically replace the zeroes to become interpolated samples?
If I decide to use polynomial interpolation (Lagrange) do I need to use a low pass filter when downsampling to remove the high frequency signals that are above the Nyquist frequency? What about during upsampling? Why or why not?
All the implementations of resampling I see on the internet are "IIR" filters. Are there any "FIR" resampling filters? If so, why is one seemingly more popular than the other?
My application works by resampling 'chunks' of audio randomly instead of a continuous stream (which is why I want to use a FIR resampling filter as there is no continuity). Is this a good idea or should I change my application to use a IIR resampling filter on a continuous stream of audio?

Thanks all

There are some excellent papers on this topic you might want to read.

Sampling-50 years after Shannon, M Unser - Proceedings of the IEEE, 2000

Splines: A perfect fit for signal and image processing, M Unser - IEEE Signal processing magazine, 1999

Interpolation and decimation of digital signals—A tutorial review, RE Crochiere, LR Rabiner - Proceedings of the IEEE, 1981.


How can I shape my time channel response, in order to have Gaussian shaped Doppler Spread? I know how to generate a Rayleigh fading channel but I have no clue how I can shape the Doppler spectrum to be Gaussian.

Any help will be really appreciated.

ps: I am coding on matlab

See  section 2.4.3 of this reference https://web.stanford.edu/~dntse/Chapters_PDF/Fundamentals_Wireless_Communication_chapter2.pdf

If the doppler spectrum has to be gaussian, the auto-correlation of tap gains should be gaussian (which is correctly mentioned in other answer but I somehow felt more details were not captured). For a coherence time $T_c$, proportional to doppler spread $1/D_s$, will be expressed as smallest amount of $n$ for which your tap at $l$ will be significantly different from tap at $l+n$. So for example, if $n=5$, gain of taps from $0$ till $4$ will be having a guassian distribution, while the gain at tap $0$ will be uncorrelated with tap at $5$.

If you model your taps as Rayleigh fading channel, your PSD will not be Gaussian as given in eq 2.60 of the reference. You may need to generate your tap gains such that their auto-correlation function follows gaussian distribution.

Since you are using MATLAB

raychan = comm.RayleighChannel('DopplerSpectrum',{doppler('Gaussian')});

This makes sure that each of the channel taps gain auto-correlation (and PSD) follows Gaussian distribution.


How to produce color image using edge detected images? Hi I am trying to produce an image which is a combination of multiple images such as in focus stacking. Now I have 6 photos and I want to produce a final image such that the final image uses blocks of individual images with maximum value of the block.
But I want to do edge detection of these 6 images first and based on the maximum edges produce a final image such that the final image has the maximum edges from the combination of all images. But the final image produced should be a color image and not in the edge detected form. Is it possible to do so?
Images can be accessed here

First of all, thanks for sharing this problem. I hope it can foster some good discussions here at stack-exchange. I guess this problem is commonly called Multi-Focus Multi-Image Fusion (MF-MIF).

Let's look at all the images together with their grayscale histograms and the gradient magnitude (information content in each image):

As the original question requires an edge-based (or rather gradient as I take it) method, I will follow the idea presented in the following paper:

  S. Paul, I.S. Sevcenco, P. Agathoklis, "Multi-exposure and Multi-focus
  image fusion in gradient domain", Journal of Circuits, Systems, and
  Computers, Vol. 25, No. 10 (2016), DOI: 10.1142/S0218126616501231
  https://pdfs.semanticscholar.org/114b/835c5ea86e93efff1027cd6da7188c0cc92f.pdf

Essentially, the method aims to blend the gradients of the luminance components of the
input images using the maximum gradient magnitude at each pixel location and then to obtain the fused luminance using a Haar wavelet-based image reconstruction technique. The algorithm also includes a Poisson solver at each resolution to eliminate artifacts that may appear due to the non-conservative nature of the resulting gradient.

The official source code of the method is given in this github repository and in your, kind of challenging dataset, the output looks like the following:

I am not sure if this is a sufficiently good output, but I bet further investigations can quickly lead to improved results. A small observation: The first image seems to be quite sharp. It seems like it has a global focus and already looks (only looks though) better than the resulting image.


How can I detect buttons in a GUI of an app with pattern detection using python? I will explain taking as example the GUI of the Calculator app in windows.
Given this Calculator application open and focused, I need to find a way of detecting all the buttons.I can only use non-intrusive methods, so things like button id are out of the question. This leads me to image recognition. Or better say image detection, because I don't want a specific image, but a set of images that follow a certain pattern. I know how I can click/right-click/dbl-click/etc using a given button image and image recognition [1]. I don't know how I can do the other way round: scan the GUI and find those areas that meet the requirements of being a button (rectangle, text/icon/graphic labeled etc.). A bigger tackle would be to detect items that don't have the rectangle shape (e.g. icons on a windows desktop)

The nearest thing to what I need is detecting faces in a picture.[2] But I don't know how to apply this in my case. For human face detection, I saw that hundreds of pictures of the face must be used in the Haar cascade generation (don't know how I would do this with only 10-15 button snapshots). If another kind of object has to be detected like an apple, you need to generate Haar cascade for that object again using a lot of images.

Did any of you ever tried to detect buttons, items or whatever in a GUI using only pattern detection? I just need something to tell me "this is an icon/button" so I can grab that region in a snapshot.

[1] I use SikuliX with Python to perform actions on given patterns.

[2] I saw that this is done easily using OpenCV and Haar cascades (in XML format). Creating the Haar cacades requires a little bit of patience and skills.

First, have a look at the squares.py sample provided by OpenCV. It should handle a fair number of button types with some tweaking.

Here is the output I got (with some tweaking) for your Calculator example:

I made the following tweaks to the squares application:

Change this code (starting on line 84):

if(result.total == 4 and
   abs(cv.ContourArea(result)) > 1000 and
   cv.CheckContourConvexity(result)):

To this:

if(result.total == 4 and
   abs(cv.ContourArea(result)) < 1300 and
   abs(cv.ContourArea(result)) > 300 and
   cv.CheckContourConvexity(result)):

Also, since you don't have to worry about scale or rotation variations, check out using matchTemplate. Also, its corresponding tutorial (in C++) is here.

See my other answer for another example of how template matching works. Also, you might find this answer helpful for detecting the top X matches with matchTemplate.

Hope that helps!


Channel noise - A misnomer? I was discussing with a colleague with background in communication about white Gaussian (WG) channel noise.

He said noise is not added by the channel and it's just a misnomer. Noise, which is thermal noise, is actually added by receiver. Channel will just introduce fading, not noise. I was stunned, as I studied in many textbooks that channel add WG noise to the transmitted signal. Is this correct?

So, it's generally right that a receiver adds noise due to physics (Johnson-Nyquist noise). So, yes, in any system, at least part of the noise in the received signal is caused by the receiver.

There's also noise that is background noise from e.g. cosmic sources that happen neither in the receiver nor the transmitter (and also aren't other transmitters interfering), but those sources are typically far less important than receiver noise – at least in the microwave bands, where modern comms happen; in LF/MF/HF bands, background noise is very observable.