What is an arithmetic mean filter? Defined as:

  $$\hat{f}(x,y)=\frac{1}{mn}\sum_{(s,t)\in S_{xy}}g(s,t) $$

As far as I understand:

$g(x, y)$ is the original image
$g(s, t)$ is a sub-image of $g(x, y)$ with a dimension of mxn
$\hat f(x, y)$ is the filtered image
Sub-images are summed up and then multiplied by ${ \frac {1}{mn}}$.

Am I correct?

If yes, then, I have several questions:

Is $g(s,t)$ a pixel in the neighborhood? What does it mean by that?
What is $S_{xy}$ ?
What are the meanings of $s$ and $t$?
What is the meaning of $(s,t) \in S_{xy}$?

In one picture: here, the mask is of size $m=3\times n=5$, filled with value $\frac{1}{3\times 5}$, the light gray shadow on the input image defines the neighborhood of the dark gray pixel: a set of pixels (all your $g(s,t)$) close to the center one (it moves with the pixel). Let us suppose the dark gray pixel has coordinates $x=5$ and $y=6$. The shadow defines $S_{5,6}$ having $s\in[5-(m-1),\dots,5+(m-1)]$ and $t\in[6-(n-1),\dots,6+(n-1)]$. $(s,t)\in S_{5,6} $ means for instance that the light gray pixel on the bottom left of the dark gray one, with coordinates $(6,5)$ belongs to the shadow.


Tricks for computing rapidly changing biquad coefficients I have 5 'negative peaking' biquads that will change center frequency while maintaining a constant Q and A (amplitude reduction) on a realtime microcontroller (ARM cortex-m4) with a codec over i2s controlled by i2c system operating at 96kHz. I am concerned (with concerns raised by my oscilloscope) that I am going to be hosed by divisions and sine/cosine. What can I do to speed things up?

Using RBJ's EQ cookbook, a peak filter takes the following coefficients:

w0 = 2*pi*f0/Fs;
A  = sqrt(10^(dBgain/40));
alpha = (sin(w0)/(2*Q));

b0 =   1 + alpha*A;
b1 =  -2*cos(w0);
b2 =   1 - alpha*A;
a0 =   1 + alpha/A;
a1 =  -2*cos(w0);
a2 =   1 - alpha/A;

Furthermore, the feedback coefficients will take the form
    c0 = b0/a0
    c1 = b1/a0
    c2 = b2/a0
and the feedforward coefficients:
    c3 = a1 / a0;
    c4 = a2 / a0;

Expanding these divisions:

c0 = (1 + (sin(w0)/(2*Q))*A)/(1 + (sin(w0)/(2*Q))/A);
c1 = -2*cos(w0) /  (1 + sin(w0)/(2*Q)/A);
c2 =  (1 - (sin(w0)/(2*Q))*A)/(1 + (sin(w0)/(2*Q))/A);
c3 = c1;
c4 = (1 - (sin(w0)/(2*Q))/A) / (1 + (sin(w0)/(2*Q))/A);

Using simplify() in MATLAB with symbolic variables w0, Q, and A, we obtain:

c0 = (2AQ + A^2*sin(w0))/(2AQ+sin(w0))
c1 = -(4AQ*cos(w0))/(sin(w0) + 2AQ)
c2 = -(A^2*sin(w0) - 2AQ)/(2AQ + sin(w0));
c3 = c1;
c4 =  -(sin(w0)- 2AQ)/(2AQ + sin(w0));

If Q and A are constant, no problem, these can be computed ahead of time. What will kill me are the divisions and the sine/cosine computations.

What tricks or processes can be used to speed up the calculation or if at all, further simplify or approximate the filter coefficients?

there is an earlier answer regarding how to more cheaply evaluate functions like sin() and log() and such.  they're approximations, but very good ones.


Why does increasing duty cycle decrease the separation between spectral components? I know that the duty cycle measures the width of the signal in percentages. But why does it decrease the signal width with decreasing duty cycle?

At 50 % duty cycle a pulse wave is a square wave. A square wave has only odd harmonics (and a possible DC bias). At other values of the duty cycle there will also be even harmonics. The emergence of even harmonics can be seen as halving of the separation between spectral components.


Finding crosses and their orientation on a given image I want to recognize table structure from some arbitrary photo image and store it in some formal notation (Let it be HTML table notation).

For example I have this blurry image as input:

After binarization I achieved this:

Now to detect structure I want to detect intersections between lines, their type (T-type, X-type, or simple corner), their orientations and their location on the image. After this I'm going to use all gained information for further joining of adjacent crosses in some structure, and translate this structure into formal representation.

The problem is in detecting crosses and their orientation:

In general these crosses may be scaled and/or rotated.
Maybe someone could help with method of solving this problem? Or maybe recommend different methods for this whole task?

That's what I have written so far:

 # -*- coding: utf-8 -*-

 import cv2
 import cv
 import numpy

 original = cv2.imread("/home/user/my_photo-1.jpg")
 grayscale = cv2.cvtColor(original, cv2.COLOR_BGR2GRAY)
 smoothed = cv2.GaussianBlur(grayscale, (5,5), 0)
 cv2.imshow("original", original)
 cv2.imshow("grayscale", grayscale)

 binarized = cv2.adaptiveThreshold(grayscale, 255, cv2.ADAPTIVE_THRESH_MEAN_C, cv2.THRESH_BINARY_INV, 7, 8)
 binarized = cv2.Canny(grayscale, 50, 200)
 cv2.imshow("binarized", binarized)
 cv2.waitKey(0)

Thanks in advance for any response/idea.

Try to use cv2.HoughLinesP (see an example) to detect the straight lines (segments). As a result you will get the segments' endpoints, so finding all the intersections and its type will be quite simple maths.
The problem will be rather with finding the proper and universal parameters.

Theoretically you could use a morphological operation called hit-and-miss transform, but this will be difficult in your case, as the lines are hand-drawn and even after thinning the crosspoints may be difficult to be properly classified (as 'T', '+' or 'corner').

I would suggest you to make some morphological preprocessing to improve the binary image you have got.


Comparing multiple $\tt armax$ models I would like to create an algorithm that, given an iddata,  estimates different armax models with different na, nb, nc, nk, and finally show the compare plot.

I know how to use the armax estimation function, but how can I store the model in a variable in a such that I can compare them later?

What's wrong with just doing sys = armax(data,[na,nb,nc,nk]); as per the Mathworks manual? sys is the variable. You could do sys[i] = ... to store several of them in an array (or perhaps sys{i}).

Then, to compare them, you could use getcov() and compare the covariances of each model.  These covariances tell you how "uncertain" the model is.

Another approach might be to look at the model residuals. If the residuals are white (uncorrelated), then the model is good.  The relatively energy / power of the residuals will be a measure of how accurate the model matches the data (the smaller the residual, the better the model).


Implementation of Filters in Python I've implemented 2 filters in Python I would like to know if they look OK or not.. I've implemented an lowpass FIR Hann Filter and an lowpass  IIR Chebyshev type 1 filter. My plots look like the following: Time amplitude and frequency spectrum:
Hann Window:

Cheby 1:

Is there something that I did wrong, that I should improve??

After roughly observing the output spectrum, I found something strange.

First, I noticed that you're coping with audio(voice) signals. Because telephones are made to work in the ranges 300 Hz to 3400 Hz [1], it seems that your IIR filter has filtered out the frequency components higher than about 1kHz. Is that what you intended to do?

Second, I noticed that your FIR filter is a kind of LPF but I cannot do any further inference, because frequency/phase responses of impulse signal have not given.
(you can create an impulse signal and then plot the spectrum as you have done.)

When it comes to the design of a filter, the spec of the filter should be discussed.

First, you need to assign the cut-off frequency (or which part of the signal do you want to filtered?). From the output spectrum of your signals, I suspect that you use two different cut-off frequencies when comparing IIR/FIR filter. If you just wanna filtered out high frequency noise, you need firstly decide the range. Maybe you can refer to [2].

Here are some frequently discussed basic spec from my experience.

Is application real-time or not?
This can determine whether you can use the so-called filtfilt technique. (Zero delay at expense of causality. (you need to recorded all signal and))
Is General linear phase (GLP) required? If so, IIR in filtfilt or GLP FIR Filter are your choices.
Usually, IIR can obtain better stop band attenuation with lower orders. However, the delay might vary from different frequencies to other frequencies. For extremely bad example, the input is 'do re mi' but the output becomes 're mi do'.

There are some more info like [3], hope it can help.

[1] https://en.wikipedia.org/wiki/Voice_frequency

[2] What is the bandwidth of human speech?

[3] https://www.gaussianwaves.com/2017/02/choosing-a-filter-fir-or-iir-understanding-the-design-perspective/


IIR filter design in matlab I have to design 2 IIR bandpass filter with the following specifications:

Sampling frequency 1000 Hz
Pass-band 50 Hz to 200 Hz
Order 6
using bilinear transform and impulse invariance with Butterworth prototype analogue filter.

I'm confused, because the two methods I've got for producing the filters have different results.

fs = 1000;
fNq = fs/2;
f1 = 50/fNq;
f2 = 200/fNq;
[z,p,k] = buttap(6);
[A,B,C,D] = zp2ss(z,p,k);
[zd,pd,kd] = bilinear(z,p,k,fs);
Bw = f2-f1;
Wo = sqrt(f1*f2);
[At,Bt,Ct,Dt] = lp2bp(A,B,C,D,Wo,Bw);
[b,a] = ss2tf(At,Bt,Ct,Dt);

[Ad,Bd,Cd,Dd] = bilinear(At,Bt,Ct,Dt,fs);
[bz, az] = ss2tf(Ad,Bd,Cd,Dd);
fvtool(bz,az);

(I think that) The above code creates an analogue prototype, transforms to a state-space form, applies the bilinear transform, transforms the low-pass filter to a band-pass filter, transforms to a transfer function and the plots it.

My problem is that it is noticably different to:

[num,den] = butter(6, [f1, f2], 's');
[B,A] = bilinear(num, den, fs);
fvtool(B, A);

I have similar problems with the impulse invariance method.

Could someone please give me a hint as to what I'm doing wrong?

The issue is resolved: the butter() function automatically does bilinear transform, you just need to specify the frequencies normalised to the Nyqyist freq.


Comparing signals recorded using diffrent methods Suppose we have two signals( focus to audio) recorded with different methods,
or more simply, suppose we have two audio signals recorded with different microphone, and  we want to confirm that these signals are of the same audio.

What are the available alogorithms or libraries that make this possible?

Different flavors of cross correlation are probably the best choice here. It really depends on how exactly the recordings are different. Two microphones side by side behave very differently than two microphones in different spots in the room. Your best method will depend on a clear understanding of what your application needs to cover and what doesn't.


Segmentation step when using feature descriptors In my final year project, we are using a SIFT feature descriptor to classify and recognize objects, but I have never come across segmentation in the process of image classification/categorization.

Isn't this step necessary?

SIFT works on points in an image while segmentation is about dividing up the image into regions. So, no, segmentation is not necessary when using SIFT.

In segmentation you divide up the image into regions so classification can be done by extracting features from each region and see if you can recognize your object. A downside to this approach is that producing a good segmentation is a difficult problem. Many segmentation techniques are sensitive to various circumstances, such as lighting, contrast, colors etc. These difficult situations often result in regions that do not represent your objects very well.

In SIFT you focus on specific points (keypoints) in the image instead of entire regions. If you want to recognize an object in an image using SIFT, you extract keypoints from an image of your object and give each keypoint a description (keypoint descriptor) based on its local neighborhood. Then you check if similar keypoints exist in the image you want to analyze. By matching the keypoints between your object image and the image you want to analyze you can do object recognition without segmentation.


Confusion about collection of bits and symbols I'm taking an intro communications module and I'm quite confused on what the term "symbol" really means?

Consider this setup:

My understanding is as folllows

After the information source, we have  a stream of bits coming in 0 1 0 1 0 0 0 1 1.....
The source encoder takes these bits and groups them into sets of N bits.
The channel encoder takes these sets of N bits and based on what coding scheme we use (gray coding) for example, it will multiply it by a certain coefficient of a orthogonal basis function (carrier signal) depending on what exact gray code it is.
The digital modulator now receives a coefficient of it's carrier/basis function and just mulitplies it by the fixed carrier/basis function (either sin or cos of something) and sends it onto the channel.

My question is - what point in this transmitter do we call "the symbol"? Is it the group of N bits after the source encoder? Or is it the coefficient of the orhtogonal basis function after the channel encoder?

In digital communications, information is transmitted in quantized form. We are constrained to transmit items that belong to a finite, discrete set. The items are called "symbols", and the set is called a "constellation". The symbols may be transmitted by changing a signal's frequency, phase and/or amplitude.

As an example, pulse amplitude modulation (PAM) is a technique where the symbols are pulse amplitudes, and the transmitted signal is a sequence of pulses. Let's say that you are using the 4-PAM constellation $\mathcal{C} = \lbrace -3, -1, 1, 3 \rbrace$, and the pulse shape is $h(t)$. Then, you group your bits 2 by 2 and assign them symbols as follows:

bits | symbol | pulse
-----------------------
 00      -3     -3h(t)
 01      -1     - h(t)
 11       1       h(t)
 10       3      3h(t)

Recall that a PAM signal can be written as $$s(t) = \sum_k c_k h(t-kT_p),$$ where $T_p$ is the pulse interval (the pulse rate is $R_p = 1/T_P$) and $c_k \in \mathcal{C}$ are the symbols. Then, transmission of the bit sequence $00111001$ would be done with the signal $$h(t) = -3h(t) + h(t-T_p) + 3h(t-2T_p) - h(t-3T_p).$$ Note that the information (bits) is conveyed by the "symbols" or pulse amplitudes (hence, PAM = pulse amplitude modulation). The task of the receiver is to recover the transmitted symbols from the received signal.

(Note that there's no unique or "official" way to map bits to symbols in 4-PAM or any other constellation; the mapping above is just an example. The only rule is to use Gray encoding, in order to minimize the bit error rate. And, or course, the mapping must be known at the receiver!)


Downsampling after a highpass filter Question: The following  sentence is not clear to me (especially why down-shift occurs). Could someone explain it please?

On the other hand, downsampling the output of the highpass filter down-shifts the high-frequency
band and expands it to the full range $|\omega| < \pi$.

This is what happens when a signal is downsampled: you get a scaling of the frequency axis and the addition of shifted spectra (aliasing). The same happens to the lowpass filtered signal. The spectrum of the downsampled signals is given by Eq. $(4.77)$ in Oppenheim and Schafer's Discrete-time Signal Processing (3rd ed), which for downsampling factor $M=2$ can be written as
$$X_d(e^{j\omega})=\frac12\Big[X\left(e^{j\omega/2}\right)+X\left(e^{j(\omega/2-\pi)}\right)\Big]\tag{1}$$
The second term on the right-hand side of $(1)$ fills up the "empty spaces" in the frequency ranges that are eliminated by the lowpass and highpass filters, respectively.


Finding horizontal count of items in an image I'm trying to find the count of horizontal columns in an ordered image (a sprite image that is used to generate an animated image):

From this image, I'm trying to generate a meta-data file like so:

{
"frames": 12,
"cols": 6,
"width": 102,
"height": 89
}

Frames would be the total count of "elements" in the image, columns is the horizontal column count, and width/height is for the individual "element" and can be calculated, once you know the column count.

What I Have Already

I've been using numpy/scipy/python mainly since I am a Python guy, but not too experienced with scipy.

This answer showed me how to get the total count of elements in an image and got me this far ipython:

Where I'm Stuck

On most simple images, the ndimage.label works well, but once the image gets too dense, it seems to fall apart.

I've tried using a kind of line-by-line average of the file: subtract 1 horizontal line, get the ndimage.label count of total elements in the picture, then keeping subtract the whole image to build buckets and calc the average std deviation for the image; which would be the horizontal count. But this still isn't accurate on dense images like this:

Would anyone have suggestions for a better algorithm or counting method?
Thanks!

For ease of hacking I have a zip with the ipython notebook and sample images available here: http://asktherelic.com/priv/sprite_sample.zip

And it looks something like this currently:
http://asktherelic.com/priv/sprite_sample.html

Using any image processing library, here is the simple step

Threshold invert  with high value will give you each object as a contour.
Now process each contour and find centre of mass.
Then sort contour by co-ordinates.


Magnitude of Power Spectral Density with Different Sampling Frequency For exponential signals (sine or cosine), if the sampling frequency $f_s$ is equal to the length of signal, $N$, the magnitude of psd for each sine signal is proportional to amplitude in time domain, while if $f_s$ is not equal to length of signal $N$, for example here, $f_s=1000 \text{Hz}$, the magnitude of sine signals in psd (frequency domain) is not proportional to their original signal. The MATLAB code is provided.

Can you tell me the reason? Is this something related with leakage?

N=256;
x0=zeros(N,1);
fs=1000;
idx = (0:(N-1))'; % Indices
f=20*[1 2 3 4];
A=[1 2 3 4];
K=4;
for ii=1:K
    x0=x0+A(ii)*exp(sqrt(-1)*2*pi*f(ii)*idx/fs);
end
X=fft(x0);

figure;
plot(abs(X).^2)

The issue (beside the scaling) here is that your selected frequencies do not align with your FFT grid and that the energy of each complex exponential is smeared out over multiple FFT bins. Your frequency resolution is given by sample rate divide by FFT length. In your case that's 3.906 Hz. Only for frequencies that are an integer multiple of this will all energy fall into a single FFT bin. An easy way to achieve this is to bump the sample rate up to 1024 Hz. Then your frequency resolution is exactly 4 Hz and all your frequencies are an integer multiple of that. Once you fix the scaling as well, it works fine like this:

N=256;
x0=zeros(N,1);
fs=1024; % Fix 1: get to exactly 4Hz frequency resolution
idx = (0:(N-1))'; % Indices
f=20*[1 2 3 4];
A=[1 2 3 4];
K=4;
for ii=1:K
    x0=x0+A(ii)*exp(sqrt(-1)*2*pi*f(ii)*idx/fs);
end
X=fft(x0);

figure;
plot(abs(X)/N); % Fix 2: absolute value only, rescaled


Extracting positive frequencies of discrete-time signal Convolution in the time domain is the same as multiplication in the frequency domain.

My data is sampled at 200 Hz, which means that the Nyquist frequency is 100 Hz, and all frequency content is <= 100 Hz. If I take a rectangle function in the time domain, with 1s everywhere below 100 Hz, FFT it, and convolve my data with that, it shouldn't make a difference. Indeed, the FFT of the rectangle function is sin(pi*x)/(pi*x), which is 0 at all integers except the origin.

When FFTing something like cos(2*pi*x/40), a sine wave at 5 Hz, you get two peaks. One at 5 Hz as expected, and one at -5 Hz. This can be thought of as two spirals e^(2*pi*i*x/40)/2 and e^(-2*pi*i*x/40)/2, one rotating counterclockwise, the other clockwise. When added, their complex parts cancel out, leaving only a real sine wave.

I want to extract only the part of my data with positive frequencies. For example, cos(2*pi*x/40) should become e^(2*pi*i*x/40)/2. This should reasonably make the data contain complex samples.

I thought to do this by creating a rectangle function with 1s where 0 <= x <= 100 Hz, FFTing that, and then convolving my data with the result. The resulting function (according to WolframAlpha) is 0.5 * e^(-pi*i*x) * sin(pi*x/2)/(pi*x/2) (graph, real part blue, imaginary part red). The problem is that the imaginary part is always zero at integers, which is weird. Since my data is real, and this function is always real at integers, the result of convolution must be entirely real, which contradicts the previous paragraph.

What is going on here? And how do I actually get what I want?

The mistake lies in the modulation of the low pass impulse response. The low pass response equals $1$ for $\omega\in(-\pi/2,\pi/2)$ and is zero otherwise, and you want a frequency response that is zero for $\omega\in(-\pi,0)$ and $1$ for $\omega\in(0,\pi)$. Consequently, you need to shift the spectrum by $\pi/2$, and not by $\pi$ as you did. The corresponding impulse response of the system becomes

$$h[n]=\begin{cases}\displaystyle\frac12,&n=0\\\displaystyle\frac12 e^{jn\pi/2}\frac{\sin(n\pi/2)}{n\pi/2},&n\neq 0\end{cases}\tag{1}$$

which can be rewritten as

$$h[n]=\begin{cases}\displaystyle\frac12,&n=0\\\displaystyle0,&n\textrm{ even},n\neq 0\\\displaystyle\frac{j}{n\pi},&n\textrm{ odd}\end{cases}\tag{2}$$

So only the sample at $n=0$ is real-valued, all samples at even indices (other than zero) vanish, and all samples at odd indices are purely imaginary.

Consequently, the real part of the impulse response is just a scaled unit impulse, and its imaginary part is a scaled discrete-time Hilbert transformer:

$$h[n]=\frac12\big(\delta[n]+jg[n]\big)\tag{3}$$

with $g[n]$ the impulse response of an ideal discrete-time Hilbert transformer.


Meaning of the Transpose of Convolution I have an $n \times n$ asymmetric convolution kernel, $k(t_1,t_2)$. $k$ is zero everywhere except for in small regions near the corners.

I also have an $n \times n$ image, $f$.

Let $*$ denote circular convolution so that the convolution operator, $K(f) = k*f$, returns an $n \times n$ image.

What is a formula for the transpose of $K$? That is, the operator $K^T$ such that $\langle K(f) ,g \rangle = \langle f, K^T(g) \rangle$.

I worked it out to be $K^T(f) = k(-t_1,-t_2)*f$ but that is not holding up to my checks.

%test how to make convolution transpose

clear all;close all;clc;

n = 1024;

%my kernel, it is complex.
[t1 t2] = meshgrid(1:n,1:n);
k = fft2(exp(-t1/50).*exp(-t2/75));

%convolution functions
convk = @(u) ifftn(fftn(k).*fftn(u));
rk = rot90(k);
convkt = @(u) ifftn(fftn(rk).*fftn(u));

%my images
f = rand(n);
g = rand(n);

trapz(trapz(convk(f).*g)) %<Kf,g>
trapz(trapz(f.*convkt(g))) %<f,Kg>, not the same as the above line...

edit:
$$
\langle k*f,g \rangle = \int_{x,y}g(x,y) \int_{t_1,t_2} k(x - t_1, y-t_2)f(t_1,t_2) dt_1dt_2dxdy
$$
$$
= \int_{t_1,t_2}f(t_1,t_2) \int_{x,y} k(x - t_1, y-t_2)g(x,y) dxdydt_1t_2
$$
$$
= \langle f, k(-1)*g \rangle
$$

edit:
Ok I got it now. There was an off-by-one error. You can see this by examining the (square) circular convolution matrix. This 1D example explains how to do it.

n = 256;

%convolution kernel
k = randn(n,1) + 1j*rand(n,1);

%circular convolution matrix
K = zeros(n);
for j=1:n
    K(:,j) = circshift(k,j-1);
end

%circular convolution fft-based function
convk = @(u) ifft(fft(k).*fft(u));

%the transpose of circular convolution, time-reversed complex-conjugated
%shifted-by-one kernel.
kt = circshift(flipud(k),1);
convkt = @(u) ifft(conj(fft(k)).*fft(u));

%test that they are the same thing
x = randi(100,[n 1]) + 1j*randn(n,1);
y = randn(n,1) + 1j*rand(n,1);
fprintf('The right matrix for K? : %d \n', norm(K*x - convk(x)));
fprintf('The right matrix for Kt? : %d \n', norm(K'*x - convkt(x)));
fprintf('The dot product test : %d \n', dot(convk(x),y) - dot(x,convkt(y)));
fprintf('The dot product test : %d \n', dot(x,convk(y)) - dot(convkt(x),y));

Will output:

The right matrix for K? : 5.246206e-11
The right matrix for Kt? : 3.582131e-11
The dot product test : 3.492460e-10
The dot product test : 9.313226e-10

edit: For 2D, I think I've got it right but I'm seeing errors around $10^{-8}$. I want to just blame round-off since if anything changes in the below code I get gigantic errors.

clear all;close all;clc;

%nxn array
n = 256;

%make convolution kernel
[t1 t2] = meshgrid(1:n,1:n);
k = fft2(exp(-t1/50).*exp(-t2/75));

%circular convolution fft-based function
convk = @(u) ifftn(fftn(k).*fftn(u));

%the transpose of circular convolution, time-reversed complex-conjugated
%shifted-by-one kernel.
kt = rot90(k,2);
kt = conj(kt);
kt = circshift(kt,[1 1]);
convkt = @(u) ifftn(fftn(kt).*fftn(u));

%test that they are the same thing
x = rand(size(k)) + 1j*randn(size(k));
y = randn(size(k)) + 1j*rand(size(k));

%check the output error
fprintf('The dot product test : %d \n', sum(sum(conj(convk(x)).*y)) - sum(sum(conj(x).*convkt(y))) );

edit: Oops, I should be looking at relative error.

abs(sum(sum(conj(convk(x)).*y)) - sum(sum(conj(x).*convkt(y))))/sum(sum(conj(convk(x)).*y))

is close to machine precision.

The transpose (or more generally, the Hermitian or conjugate transpose) of a filter is simply the matched filter. This is found by time reversing the kernel and taking the conjugate of all the values.

Chapter 5 of Jon F. Claerbout - Earth Soundings Analysis: Processing versus Inversion has a great discussion of the topic of Adjoint operators (or transpose, or back projection, or matched filtering). It's in the context of Seismic imaging, but it really gets to the heart of what you're actually trying to do (in general) with these techniques.


Magnification of a colored image through two color space: RGB & YCbCr I have tried to magnify a colored image using bicubic method. I did this in two ways; First, I resize R, G, and B parts of it. Second, I transformed RGB to YCbCr color space and then resized Y, Cb, and Cr parts of it. When I use YCbCr color space, final image changes in color but not in whole shape of the image. I could not realize the reason. May you help me?

code:

clear all
clc

A=imread('Sun.jpg');
R=A(:,:,1);
G=A(:,:,2);
B=A(:,:,3);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Y  = 0.257.*R + 0.504.*G + 0.098.*B + 16;
Cb = -0.148.*R - 0.291.*G + 0.439.*B + 128;
Cr = 0.439.*R - 0.368.*G - 0.071.*B + 128;

Y1  = imresize(Y, 2, 'bicubic');
Cb1 = imresize(Cb, 2, 'bicubic');
Cr1 = imresize(Cr, 2, 'bicubic');

R1 = 1.164.*(Y1-16) + 1.596.*(Cr1-128);
G1 = 1.164.*(Y1-16) - 0.813.*(Cr1-128) - 0.391.*(Cb1-128);
B1 = 1.164.*(Y1-16) + 2.018.*(Cb1-128);

A1 = cat(3,R1,G1,B1);

imwrite(A1,'Sun1.jpg');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
R2=imresize(R, 2, 'bicubic');
G2=imresize(G, 2, 'bicubic');
B2=imresize(B, 2, 'bicubic');

A2 = cat(3,R2,G2,B2);

imwrite(A2,'Sun2.jpg');

Original image:

Resized image(YCbCr color space):

I have checked your code and it worked fine in my matab 6.0 implementation with the following lines of of modification:

A=double(imread('Sun.jpg'));  % convert uint8 data to double
...                           % further operations here
A1(A1(:)<0.0) = 0.0;          % check and correct out of range values
A1(A1(:)>255.0) = 255.0;      % check and correct out of range values
imwrite(A1/255.0,'Sun1.jpg'); % while writing the data of class double use [0,1] data range

According to this, as I said you must deal with numerical issues such as an out of range value.


LPF - signal values unaffected at specific times Is it possible to design an LPF that has an output identical to the input at specific points in time domain (the rest of the input waveform can get filtered/distorted)?
Is there a general name/technique for this kind of thing (assuming it is possible), so that I can search for more information on this topic?

What you probably are looking for are Nyquist/$M$th-band filters. The impulse response of these have the following property (assuming a non-causal impulse response centered around tap 0 for ease of exposition)

$h_n = \left\{ \begin{matrix} \frac{1}{M} & n=0\\ 0 & n = kM, k=\pm 1, \pm 2, \dots\end{matrix} \right.$

This means that every $M$th sample is untouched (the $1/M$ term is for scaling purposes). Often it is said that these filters have zero inter symbol interference, assuming that a symbol arrived every $M$th sample.


Algorithm to extract polygones with holes I have this picture and would like to extract the bones.
First I thought of contour-Tracking, but as you can see in some bones there are holes (Like for example at the left finger, in the bone in the middle).

What algorithm would you recommend me? To find the contour of the bones? How could I fill them? Thanks

Im also grateful for code examples.
.

You need to apply some image segmentation algorithms. You may read about A Comparison of Image Segmentation Algorithms and apply the ones that satisfies your needs better.

If you're developing with OpenCV, you can find some ready to use (implemented) algorithms that can save you lot of work as it is the case with the watershed algorithm.

Note that you can use OpenCV with the traditional programming language C/C++ but also many others such as Python and Java. Personally I experienced it with C/C++ and Python though.


How to accurately model wideband RF reflections? I would like to model a transmitter and received pair, both sampling at the same bandwidth (not clock synchronized), where there are a number of reflectors. The time delay due to the reflectors can be a non-integer multiple of the receiver resolution (1 / sampling bandwidth), and I would like to accurately model that.

This means that if I take the transmitted signal at the nominal sample rate, and sum a bunch of time-delayed versions of the signal (with the time delay equal to the difference in time it takes light to travel the direct path and the reflected path), this will be inaccurate because the time delays will be quantized to the sampling period, which is too coarse.

Obviously I can upsample and then do the sum of the upsampled discrete-time signals, which will decrease the quantization error of the time delay. But I'm wondering if there's a way to model the received digital signal without any quantization error. In other words, I'd like to model the continuous-time reflections.

I'm modeling this for a RADAR application, so I'd like the model to be phase-accurate. I'm assuming that each reflector can be modeled by the distance to the transmitter and the receiver, and some real-valued attenuation constant. Is this assumption accurate for RADAR applications? I'm also assuming that the received signal is a sum of time-delayed signals with time-delay and attenuation from the reflector's distances and attenuation constant. Is this accurate for RF, or do I need to account for a discrete phase jump at each reflection?

So, to give you something to read up on first, the channel you describe is a Rician or Rayleigh channel, depending on whether you have a dominant line-of-sight path or not.

So, as a first approach, to delay something in time, you don't have to shift it by a whole sample – you can also do it in frequency domain, by DFT'ing your signal, multiplying it with a $e^{j\frac{\Delta t}{f_s} f}$, and then IDFT it back. That's how I've built all my radar simulations.

  I'm assuming that each reflector can be modeled by the distance to the transmitter and the receiver, and some real-valued attenuation constant.

A radar target typically also phase-shifts the electromagnetic wave – on perpendicularly hitting a metal surface, by 180°, but differently for more complex surfaces. In reality, we often model that stochastically. What is appropriate for you will very much depend on what you simulate. A model for scattering of reflections e.g. on leaves will not be very useful for modeling reflections e.g. inside an oil tank, and vice versa. Both are use cases for radar.

For many scenarios, there's "standard" models, such as the COST-Hata urban model, the ITU indoor models and many, many more. It really depends on which environment you're considering, and also, which group of researchers you ask.

Many of these models are defined for what I call the "easy RF" ranges below ~5.8 GHz. Since radar, for Doppler and wavelength reasons, often operates above, things often aren't too well represented by these models (i.e. they simply don't apply to common radar frequencies > 8GHz). This has two interesting aspects:

The higher your frequency, the more quasi-optical your models get for most useful spatial scales, and
Since millimeter waves are very en vogue right now (with all the 5G research going on), you might find very interesting RF propagation models when looking at very recent papers. As an example for a very practical, hands-on-testing based propagation/reflection description, try Simić, Perpinias and Petrova: 60 GHz Outdoor Urban Measurement Study
of the Feasibility of Multi-Gbps
mm-Wave Cellular Networks. Downside of these considerations is always that, albeit multipath does play an important role for these researchers, the individual phase really doesn't matter (only the interference with other paths does).


State space equation from differential equation I have very general system. I don't know whether it is electrical or mechanical or whatever. This system can be modeled by the following differential equation

$$\dot q = \frac{Tf_1-f_2}{T+1}$$

where:

$\dot q$ would be equivalent to current in an electrical system
$f_1$ and $f_2$ would be equivalent to current sources
$T$ is a given constant

From what I see, if $f_1$ and $f_2$ are constant, there is no way for $\dot q$ to change. How is it possible that this system has a state space representation? And how do I get this state space representation?

Assuming that the inputs are $f_1$ and $f_2$,

$$\dot q = 0 \cdot q + \begin{bmatrix} \quad\left(\frac{T}{T+1}\right)\\ - \left(\frac{1}{T+1}\right)\end{bmatrix}^\top \begin{bmatrix} f_1\\ f_2\end{bmatrix}$$

Since you have not specified the outputs, the output equation cannot be written.


Comparing arithmetic complexity of FFT radix-2 and convolution Let's assume we have a discrete linear time invariant system and we have a real signal $x[n]$ with length N=50 as input for the system.

The impulse response $h[n]$ of the system is considered to be also a real signal with length M=10.

We want to calculate the output signal $y[n]$.

There are two ways we can do this: We can either convolute the signals $x[n]$,$h[n]$ in the time domain or use FFT radix-2 and work on the frequency domain.

I want to compare the arithmetic complexity of these two methods in number of real multiplications.I know that the arithmetic complexity of convolution is $N \cdot M$ since $N \cdot M$ real multiplications are required in order to calculate $y[n]$. However, I am bit confused with FFT-radix 2.

I know that it requires $μ(Ν) = \frac{N}{2}\log _{2}N$ complex multiplications. Every complex multiplication is equivelant to 4 real multiplications. So how many real multiplications are required to calculate $y[n]$? Do I have to calculate first both $X[k]$ and $H[k]$ by using $\frac{50}{2}\log _{2}50 + \frac{10}{2}\log _{2}10 $ real multiplications? And if so, how do I go on?

Let the input signal $x[n]$ has a legth of $N=50$ samples, and the filter $h[n]$ has a length of $M=10$ samples, then the output $y[n]$ (by linear conv) will have a length of $L = N+M-1=59$ samples.
If you use time-domain convolution, then the number of real MACS can be seen to be  $N \cdot M = 50 \times 10 = 500$.
If you want to use radix-2 FFT to implement the linear convolution result, then you should select a length of $R = 64$ for FFTs. And you will: 1- convert $x[n]$ and $h[n]$ into $X[k]$ and $H[k]$ by two $R$-point FFTS, 2- multiply the results to get $Y[k] = X[k]H[k]$, and 3- apply inverse FFT of $R$-point on $Y[k]$ to get the output $y[n]$.
Each $R$-point FFT (and IFFT) requires about $\frac{1}{2} R  \cdot \log_2(R)$ complex MACs. One complex MAC is 4 real MACs, therefore this is equivalent to $ 2 ~R ~ \log_2(R) = 128 \log_2(64) = 128 \times 6 = 768 $ real MACs. And the total of two FFTs and one inverse FFT requires about $3 \times 768 = 2304$ real MACs. The intermediate multiplication also requires $4 \times R = 256$ real MACs and hence FFT based implementation requires a total of $2560$ real MACs.
So in this case, time-domain convolution is more efficient than an FFT based implementation.


3dB and 6dB Per Octave Rolloff Digital Filters / Generating Pink and Brown Noise So I would like to have a design method for filters with 3dB and 6dB per octave roll off -- for generating pink and brown noise respectively.

I know the following 'pinking' filters exist:
Filter to add 3dB per octave?

But the poster says he doesn't remember the value for A, nor does he give any idea of how he did it in the first place.

So my main question is:

How do we pick the poles and zeros of such a filter?

-6 dB / octave is trivial.

a common reference for a pinking filter is http://www.firstpr.com.au/dsp/pink-noise/


Mel Cepstral Distortion I am working on a speech synthesis model and I am looking to evaluate my synthesized speech. I found that most people use the Mel Cepstral Distortion (MCD) which can be calculated by the formula:

$$\frac{10\sqrt{2}}{\ln10}\frac{1}{T}\sum_t\sqrt{\sum^{25}_i 2||mc(t,i)-mc_{synth}(t,i)||^2}$$

I have ground-truth speech and the synthesised version so I calculate the MFCCs for both (taking 25 coefficients) and then apply the formula. The result I get is around 300, which seems orders of magnitude larger than what I should get. I tried doing cepstral mean normalization which brings the number down to around 30 but I notice that most of the time this number is <10. The speech I synthesize is not so bad so I must be doing something wrong. Does anyone have any ideas?

P.S. The way I produce the speech ensures that it be aligned with the original(more or less). I wonder if a slight misalignment could cause this issue and if so should I use DTW to align the audios?

Edit
Also If anyone knows of a library that calculated MCD when given an audio file then I could use that to compare my results.

Ok so here is what I found. The distance is dependent on the way that the mfccs are calculated. This makes sense to me and also explains why cepstral mean normalization affects the values of the MCD.

I found this implementation (https://github.com/MattShannon/mcd), which unfortunately did not support .wav files. I ran this and it gives results that are in the correct range for its test examples. I looked into the code and found that the the difference was not so much in the calculation of the MCD (which is straight forward) but rather the mfccs. The library I was using gave mfccs at a different scale (I tried python_speech_features and pyspeech). Using the mfccs from Matlab gave me mfccs in a similar range to that of the repo so I opted to use those.

A noteworthy point is that the first coefficient is quite large which is why it is ignored in (https://github.com/MattShannon/mcd) so I did the same. There are some artifacts in the speech I synthesise which likely accounts for why my MCD is still above 10 but its still lower than 30 now.


What is hidden in DFT below frequency resolution? Setting: signal sampling rate 1000Hz, DFT 1s 1000 samples. This yields frequency resolution 1Hz.
I know, this is the limit of meaningful frequency and when I want to read 0.5Hz then either sampling rate should be doubled or samples acquisition doubled.
What bothers me is that when I compute DFT of 0.5Hz in given setting on artificial signal (sparse signal, some k Hz cosines and some 0.5 Hz cosines, e.g. 2Hz, 5.5Hz, 12Hz, 14.5Hz etc.) without any spectral window, the results seems valid.
I know this shouldn't happen, with dense signal, non-perfect frequencies (multiples of frequency resolution) and noise situation gets complicated.
What exactly I read below frequency resolution and what are conditions it may work?

Computing your DFT with a 0.5 Hz spacing is equivalent to appending your signal with 1000 zeros and then computing a 2000 point DFT. It should be clear that this adds no real information to the signal. Yet it can sometimes be useful.
What really happens by this zero padding is that you're seeing an interpolated version of your spectrum. This means that between each two "valid" points (at 1 Hz spacing) you're seeing one more that produces an interpolate, a potentially smoother version of your spectrum. The interpolation kernel is essentially a sinc function (periodified sinc to be accurate).
Depending on what you want to do with the spectrum, such interpolated results can be visually more pleasing and sometimes more clearly show things we would otherwise have missed.
In particular, if you have one harmonic signal and its frequency is not on the DFT grid, your DFT will not show one line but many. This is also known as the picket fence effect. This picture may be confusing and in cases like this, an interpolated version would show us more clearly what is going on and thus be easier to interpret. Also, this allows us to find the actual frequency with a better precision (NB: resolution is not improved by this, but precision is).
Here are two examples. In each example I considered N = 30 for clarity. For the first example we have a frequency of 5 Hz, falling exactly on the DFT grid. The 30-point DFT (depicted in red) gives just one line). An interpolated version is added with the dashed line just for clarity. You can clearly see the sinc-like interpolation kernels and how in the on-grid case we're sampling exactly the zero crossings, which is the reason we're seeing only one line.

In the second example, the frequency was chosen to be off the grid (between 5 and 6 Hz). Now the 30-point DFT shows many lines and it may be hard to see that these stem from just one harmonic if we only have those. Adding the interpolated DFT once more shows a very similar picture to before: We can see that it looks like one harmonic and we can much better spot its maximum.


Using DFT to determine control system bandwidth from user input I am currently working on a race car project where steering and braking is to be done by a fully electronic x-by-wire system. To roughly determine the minimum bandwidth of the controllers, I took user + car data from an existing race car. I though I would do an FFT of the input signals to see what the 'user' bandwidth is, and let the controller sample time be at least 10 times this value.

As an example, this is what I did with the steering input:

I took the sampled time-series of the steering wheel angle, s(t), for 2 laps, and removed the mean value from the series.
I determined |s(t)|^2 to obtain a power-like quantity from the user input
An FFT was performed on the signal from step 2, and the bandwidth of this signal was determined. At first I was going to determine the half power bandwidth, but I quickly noticed the spectral components above 2 Hz were negligible so I took this as the 'user signal bandwidth'
From step 3, the min sample frequency for the steering input is 20 Hz as rule of thumb.

My questions are:

was I correct in thinking that an FFT of the input signal represents
how much input 'energy' is present at each frequency in the FFT?

does step 2 makes sense or should I skip it?
Am I correct in thinking the system will steer nicely if I sample the input signal at at least the frequency from step 4?
My grad prof asked me 'so an FFT of the steering input, what am I seeing here? What does this graph actually mean? Is the driver
actually turning the wheel at 1/3 Hz, what have you plotted?' and I
had no idea how to clearly answer this question. Any thoughts?

I am new to this forum so hope my questions are clear enough for you guys.

Thanks in advance for helping me out!

Your basic approach seems correct! However, step 2 is not a good idea since taking the square of the signal is a non-linear operation. Taking the square of the FFT to determine the signal power is correct though.

For your profs. question it depends on what you have measured and at which sample rate. I assume it is the steering wheel angle and therefore it should be possible to answer it, but there will be some pitfalls as most likely the steering movement is not in a sinusoidal form. However, the basic frequencies in the movement will be related to this.

Given that you redo the estimation without the squaring of the signal, the ten times safety margin should be more than enough.

To further convince you and your prof. I would suggest bandlimiting the measured data to your estimate by using a very sharp FIR filter and compare visually as well as downsampling the data and compare.


Compressive Sensing and Sparsity We apply compressive sensing to reconstruct a signal if it is sparse in the original domain or has a sparse represetation in some basis.

How we may know a if a signal is sparse or has a sparse represetation ?
Is a radar received signal sparse or has a sparse representation ?

It depends on the application and on what basis you decide to look at. If there are only a few targets and little to no clutter, then a radar image can be considered sparse in the image domain i.e. it has only a few Fourier components. A possible example of this is ground to air or air to air radars.

Synthetic Aperture Radar (SAR) imaging of the ground - imaging point-like targets on a specular reflecting surfaces (road, highways, airport runways) may also be considered sparse in the image domain. Once you have speckle though (e.g. high res imaging of fields, trees etc) then the image is not sparse in the target domain and sparse reconstruction will have a hard time with this type of data. The spare reconstruction may be able to separate point like targets from the clutter because it reconstructs the targets and not the speckle background.

Sparse reconstruction has also been applied to Inverse-SAR imaging of ships at sea. The ship tends to behave as a collection of corner or other simple reflectors - the reflection off the ocean waves tends to be away from the radar. You do run into cases where the ocean surface is rough and you get significant returns from it.

Yet another approach in SAR imaging to assume that the gradient of the image is sparse. This tends to behave like an edge detector.

So, it really depends on the application of the radar and on what basis you decide to use.

The CoSeRa workshop is held almost every year and focuses on compressive sensing. Many of the papers are focused on radar applications (SAR, moving target indication, tomography). The papers from the last few years are available through the IEEE.


Why are zero values added in the FFT of a concatenated noise signal? After concatenating an 11-second *.WAV (sample rate 44100 Hz) file 27 times (just by gluing them head-to-tail), I obtained a ~300-sec noise signal. I used MATLAB(2015b) and the FFT function on both the original file and the concatenated signal. The FFT of the concatenated signal was, as expected, pretty similar to the original. However, for every y-axis point (spectral power) in the concatenated FFT, there were 2 zeros (not precisely zero, but values in the order of x*10-19, so about 18 orders of magnitude below the 'real' FFT values). The original noise file was specifically designed to allow for looping. Below a close up of the concatenated FFT showing the added zeros:

Hence, the spectral power plot of the concatenated signal looked like the original, but every frequency power point was preceded and followed by two zeros. How can this be? Does it have to do with the fact that lengthening a signal increases the resolution of the FFT and that frequencies are analyzed in the concatenated signal that are not present in the original signal due to resolution restrictions? Or is this rather an error in my script? Other spectra showed up fine, though.

You're running into a property of the DFT that is usually used in the opposite direction: stuffing zeros between samples in one domain results in replication of the entire sequence in the opposite domain. Let's start in the frequency domain with the signal you plotted, $X[k]$, and look at the inverse transform to see what time-domain signal it corresponds to.

$$
x[n] = \sum_{k=0}^{N-1}X[k] e^{j 2 \pi k n/N}
$$

$N$ is the length of the DFT that you apply to the data. Let's say that $X[k]$ consists of sparse values with $M$ zeros between each one, so that $X[k] = 0\ \forall\ k \neq kM$. The above then looks like:

$$
x[n] = \sum_{k=0}^{\left \lfloor{\frac{N-1}{M}}\right\rfloor}X[k] e^{j 2 \pi k M n/N}
$$

Do some rearranging within the exponent to get:

$$
x[n] = \sum_{k=0}^{\left \lfloor{\frac{N-1}{M}}\right\rfloor}X[k] e^{j 2 \pi \frac{k n}{\frac{N}{M}}}
$$

This is nothing but a $\frac{N}{M}$-point inverse DFT of the nonzero samples in your frequency-domain signal $X[k]$. Due to the periodicity property of the DFT, you can then show that $x[n]$ is periodic with period $\frac{N}{M}$, which is exactly the form of signal that you shoved into the DFT to begin with.

Summary: Replication of the signal in one domain results in a zero-stuffed signal in the opposite domain. This is used in the theory behind interpolating filters and other multirate processing. There is a dual property where if you downsample a signal (by removing $M$ out of every $N$ samples), the data in the opposite domain collapses on top of one another and adds (resulting in aliasing).


Are There Common Values of Standard Deviation for Gaussian Noise of an Image? For example, it occurs to me that there might be a standard deviation value of Gaussian noise in which the human eye could no longer distinguish an object because the standard deviation is very high and many edges were removed because of that amount of noise. Or I think maybe there are values of standard deviation in which the human eye can better distinguish an object in an image because the model of vision of the eye is modeled in such a way that these values of standard deviation allow the eye to recognize or distinguish an object. Is there any research about this that I am commenting on? Are there "important" values of Gaussian noise standard deviation in an image?
Thanks in advance

You can easily have a look on the values of the STD on Image Denoising Papers:

The range of 1-15 is considered low.
The range 15-30 is considered medium.
The range 30-50 (Even above) is considered high.

The above values is for images in the range [0, 255].


Does "improper" imply that a system cannot be stable and causal? This answer and the comments in it made me wonder whether the following statement is true:

  If a transfer function is improper, then that system cannot be causal and stable at the same time.

I had thought that this was true for a while. But the other day I wondered why. For example, the transfer function

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

is improper, but the ROC $\{s:\mathrm{Re}(s)>-1\}$ would make it stable (it contains the imaginary axis) and causal (it consists of a left-sided plane).

So how are causality and stability related in improper systems?

An improper system cannot be causal and stable. If the order of the numerator is greater than the order of the denominator, you'll always have at least one pole at infinity. Consequently, not all poles are in the left half-plane (or inside the unit circle in the case of discrete-time systems).

The system in your example is clearly unstable:

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

Part of it is an ideal differentiator ($s$), which is unstable.

Also take a look at this related answer.


Explain Auto-Tune in a simple way I have to do a presentation about Auto-Tune and its relation to the Fourier transform. What is a good explanation on how does Auto-Tune work?

when I read the original autotune patent few steps made me understand that everything was done in time domain (in the past, I don't know today), they didn't mention anything about overlap and add, it made me wonder if the pitch detector is so good that they didn't have to overlap and add, could they always skip or add periods in the exact position? (just out of curiosity I ended up implementing the pitch detector described in this patent, and honestly that's not all I imagined lol)

Take a look in the red marked piece of the patent paper:

Of course, the patent description also says nothing about keeping the formants, but to me it is clear that the above description is just a time domain method to apply time scale and change the pitch using resample.

By the way, I wrote a pure time domain "pitch correction" keeping the formants based in the papers from @robertbristow-johnson and Keith Lent

youtube demo from my code here:

https://www.youtube.com/watch?v=6ns5K1FHtd4


How to reconstruct the signal which is subjected to FIR filtering - Matlab I have a signal $x$ which is corrupted by Gaussian noise.
The signal $x$ is convolved with a filter with coefficients $h = [1, -2, 1]$, to generate signal $y$. Now the signal is taken to Discrete Fourier Transform (DFT) domain by FFT.

In MATLAB,

y = conv(x,h);
z = fft(y);

How can I compensate for the modulation caused by $h$ in the DFT domain.
I have attenuated the frequency components of $x$ by the filter $h$.

I know that to invert an FIR filtered signal, IIR filters are required.
Can I do some compensation (multiplying coefficients of $z$) on Fourier domain.

You can if the filter does not introduce any zeros in the frequency domain, which your filter unfortunately does because its frequency response is zero for $\omega=0$. What you normally have to do is zero-pad the impulse response $h(n)$ to the same length as the result of the convolution, i.e. $y(n)$. Then you compute the DFT (FFT) of the two zero-padded sequences:

$$H(k)=DFT\{h(n)\}$$

$$Y(k)=DFT\{y(n)\}$$

Now you have

$$Y(k)=H(k)X(k)\tag{1}$$

where $X(k)$ is the DFT (FFT) of the zero-padded input signal (such that it has the same length as $y(n)$). From (1) you can compute $X(k)$ for all indices $k$ for which $H(k)\neq 0$. In your case $H(0)=0$, so you can't compute $X(0)$. You can, however, compute all other values of $X(k)$ from (1).

For your example this means that you can determine the original time signal $x(n)$ from $X(k)$ up to its mean value. E.g., if you choose $X(0)=0$ and apply an inverse DFT you'll get the original signal but with a mean value of zero.

In Matlab it would look like this:

Nx = 100; x = randn(Nx,1);    % some input signal
h = [1,-2,1]'; Nh = 3;        % impulse response
y = conv(x,h);                % convolution
h0 = [h;zeros(Nx-1,1)];       % zero-padding
Y = fft(y); H = fft(h0);
X = Y(2:end)./H(2:end);       % undo filtering except for DC
X = [0;X];                    % add arbitrary DC value
x_ = real(ifft(X));           % reconstruct input signal (up to mean)
plot(x-x_(1:Nx));             % should be constant


FIR filter with relaxed transition band Is there a way to design a lowpass FIR filter with a more relaxed transition band than what MATLAB FilterDesigner tool generates?
What I intend to do is to reduce number of taps needed to implement the filter. The response that MATLAB generates is good but it needs too much coefficients and the response that I need can be relaxed in the almost second half of the transition band.
The spec. of the desired filter is as below:
Fpass = 5 M;
Rejection at 700k offset from band edge: 15 dB
Rejection at 1.5M offset from band edge: 30 dB
Rejection at 5.5M offset from band edge: 70 dB
Here are the parameters I have used (Units are MHz and dB).
Fs = 187.5;
Fpass = 5;
Fstop = 10.5;
Apass = 0.35;
Astop = 70;
Other parameters are method = equiripple and density factor = 20;
Though the resulted filter has 70dB rejection at 10.5 M, I doesn't have the rejections needed at 700k and 1.5M offsets. To get that 15 dB rejection I have to either increase rejection at stopband or decrease Fstop which will result in more coefficients and is overkill for rest of the filter.
Also tried firpm() and using the parameters below I got the same response as the one by FiterDesigner which didn't help.
>> [n,fo,ao,w] = firpmord([5e6 10.5e6],[1 0],[0.01 0.0002],187.5e6);
>> b = firpm(n,fo,ao,w);

Instead of using the FilterDesigner Tool consider using the firls() and firpm() commands directly since they will return optimized filters (in the least squares sense or peak error /equiripple sense respectively) very easily with one MATLAB command. Type help firls or help firpm in MATLAB to read further documentation on how to use these and specifically how the parameters define the transition bandwidth desired.
This can equally be done through the FilterDesigner but I see no need for that added complexity.


Butterworth filter in python I'm trying to use a Butterworth filter in Python as described in this thread with these functions:

def butter_bandpass(lowcut, highcut, fs, order=5):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a

def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    y = lfilter(b, a, data)
    return y

The output of the FFT of my data without applying the filter gives the following plot:

However, after applying the filter above with:

lowcut = 1.0
highcut = 50.0
x2_Vtcr = butter_bandpass_filter(x_Vtcr, lowcut, highcut, fs, order=4)

where fs is the sampling frequency (250 in my case) I get as FFT:

It looks like the filter shifts the frequency to the left and I don't get the peak where it should be.
Any idea as to why this is happening?

Thank you!

(Not an answer, but can't post code in a comment)

Works fine for me:

from

 __future__ import division
from numpy.random import randn
from numpy.fft import rfft
from scipy.signal import butter, lfilter
from matplotlib.pyplot import loglog

def butter_bandpass(lowcut, highcut, fs, order=5):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a

def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    y = lfilter(b, a, data)
    return y

fs = 250
lowcut = 1.0
highcut = 50.0
x_Vtcr = randn(10000)
x2_Vtcr = butter_bandpass_filter(x_Vtcr, lowcut, highcut, fs, order=4)
loglog(abs(rfft(x_Vtcr)))
loglog(abs(rfft(x2_Vtcr)))

What do you get that's different from this?


What are your tips for DSP self-study? Perhaps a meta-question, but I'll post it here because some of you are bound to possess gold nuggets about continuing education. It might be closed for being subjective but I don't care because this community is small and this could be valuable info for visitors.

Aside from reading DSP.StackExchange questions and answers, I've been doing the following:

MIT OCW Signals and Systems, but Oppenheim's lectures from the 80's and not the 2004 course, actually doing the problem sets, reviewing
the solutions, and really understanding where a process comes from
Reviewing circuit analysis diagrams
Taking the concepts from Signals and Systems and plotting / manipulating them in MATLAB

Of course, there's no substitute for working on a project and having your *** on the line for not making it work, but things can only get better if your fundamentals are really solid, right?

Tips for DSP self-study huh. Well, ...studying 'signals and systems' is a great idea and having Matlab software means you have the tools to learn an awful lot. I think Dr. Steven Smith's book "The Scientist and Engineer's Guide to Digital Signal Processing", which you can read online for free, is a terrific source of fundamental DSP information. Dr. Smith is a good writer. His book is at: http://www.dspguide.com/pdfbook.htm. Smith is a skilled working DSP engineer. That means he knows which end of the soldering iron is hot.

With regard to a free DSP textbook homework problems, I recently learned that Sophocles J. Orfanidis, the well-known professor with the ECE Department of Rutgers University, has made his signal processing textbook available for downloading on the Internet. The textbook is: "Introduction to Signal Processing" available at: http://eceweb1.rutgers.edu/~orfanidi/intro2sp/. Happily, also available at the above web site are:
(1) Errata for the Orfanidis textbook, (2) Orfanidis' Homework Solutions Manual, and (3) the errata for Orfanidis Solutions Manual.

At the risk of appearing tacky, you might consider my "Understanding Digital Signal Processing" textbook. It's been Amazon.com's top-selling DSP book for some years now. My book's been popular not because I'm any more DSP-skilled than other book authors---I can assure you that is not the case---but rather because I have some skill at explaining signal processing concepts from a working engineer's viewpoint rather than a university professor's strictly-mathematical viewpoint. My book costs money. However, I can assure you, I use the royalties from my book to buy milk and bread for my children.

So Anthony, learning the fundamentals of DSP through reading, performing lots and lots of Matlab modeling, and maintaining your enthusiasm, to quote Susan B. Anthony, "Failure is impossible." Good Luck.


Synchronization with a specified sequence I'm currently creating a project in Matlab where I'm simulating a communication, based on the SSB modulation, between a transmitter and a receiver.

I've added a Barker sequence in the trasmitted message in order to Synchronize the receiver. The receiver starts to demodulate the message from the istant when the cross-correlation beetween the Barker sequence and the received message is maximum.

The project works with no problems, but I have a question...which kind of autocorrelation function should have a sequence used for this scope? I've found that a good sequence should have an autocorrelation function similar to a pulse...is it correct? If it is...why?

Thanks for the attention!

I've found that a good sequence should have an autocorrelation function similar to a pulse...is it correct? If it is...why?

You want to pinpoint one specific delay – the right time synchronization "point". From that, it directly follows that the ACF should be as "dirac-y" as possible to minimize ambiguities and maximize energy in the one right time.


Octave - freqz strange results I wrote octave code to find transmittance based on zeros and poles:

sps = 44100;     % sample rate [samples/s]

function [b, a] = dbtf(sps, z, p)
  b = [ 1 ];
  a = [ 1 ];
  for i = 1:size(z)(1)
    f = z(i, 1);
    g = z(i, 2);
    b = conv(b, [1, -2*g*cos(2*pi*f/sps), g^2]);
  endfor

  for i = 1:size(p)(1)
    f = p(i, 1);
    g = p(i, 2);
    a = conv(a, [1, -2*g*cos(2*pi*f/sps), g^2]);
  endfor
endfunction

z = [ 0 1; 30 1;            45.8 1;          66 1; ];
p = [           41.7 0.9999;       60 0.9999;     95.2 0.9999; ];
[b, a] = dbtf(sps, z, p);
[h,w] = freqz(b, a, 128000);
semilogx(w/pi+0.00001, abs(h));

It just calculates conjugate zeros and poles of numerator and denominator. Next zeros and poles are added by polynomial multiplying using conv function. Z and p vectors holds frequency/radius pairs.

For such defined zeros and poles I have strange result:

z = [ 0 1; 30 1;            45.8 1;          66 1; ];
p = [           41.7 0.9999;       60 0.9999;     95.2 0.9999; ];

but if I remove zero at 30Hz, the result is as expected:

z = [ 0 1;                  45.8 1;          66 1; ];
p = [           41.7 0.9999;       60 0.9999;     95.2 0.9999; ];

The question: Is it just matter of number precision in octave or my calculating method is wrong ?
Alternative question: Could somebody try it on matlab (not octave) and say if results are the same ?

I of course expect notches at frequences of zeros and resonances at frequences of poles. Only second result shows 3 resonances and 3 notches as expected. In the first example frequency response looks completly random (nothing resembles expected notches and resonances).

What you are seeing are indeed numerical artefacts. The conversion from polynomial roots to polynomial coefficients (and vice-versa) is known to be an ill-conditioned problem. In your case the problem occurs at a relatively low polynomial order because zeros and poles are all clustered around the same point (DC), and the poles and zeros almost cancel each other. By removing one zero, as you tried, the problem's numerical condition is improved, simply because the polynomial order becomes smaller. Another way to avoid the problem is to distribute the poles and zeros more evenly around the circle in the complex plane. Note that the available frequencies extend up to $f_2/2=22050\text{ Hz}.$ That means that just decreasing the sampling frequency would also improve the numerical condition.

Apart from the numerical effects, your approach for converting zeros and poles to numerator and denominator coefficients is correct. However, it can be done in a much simpler way by using the function zp2tf, which transforms given poles and zeros to transfer function coefficients.

Also note that with your approach you generate a double zero at zero frequency because you produce complex conjugate pairs, even for real-valued zeros.

Finally, the filter you computed is not causal and stable. This has nothing to do with the effects you see, but I believe it's still worth mentioning. So you wouldn't be able to implement that filter with the coefficients you computed. The reason for this is that the transfer function is not proper, i.e., its numerator order is greater than the order of its denominator. This means that apart from the given poles, there occur a number poles at infinity, and that number equals the difference between the numerator and denominator orders. This problem can be solved by adding poles at $z=0$, i.e., by multiplying the denominator by $z^K$ where $K$ is the difference between the polynomial orders of numerator and denominator. This will not change the frequency response, it will just make the filter causal and stable.

EDIT:

For maximum precision you shouldn't compute the coefficients of the numerator and denominator polynomials. Instead you should work with second-order sections. Below is a simple example in Octave/Matlab. Note that I didn't optimize the sequence of sections, neither did I optimize the pole-zero pairing per section. But, as you can see in the figure, the numerical artefacts are greatly reduced.

fs = 44100;
rz = 1; rp = 0.9999;
phi_z = 2*pi*[30; 45.8; 66] / fs;
phi_p = 2*pi*[41.7; 60; 95.2 ] / fs;

m = length(phi_p);
L = 2^15;

a = zeros(m,3); b = zeros(m,3);
H = zeros(m,L);

for k = 1:m,
    b(k,:) = [1, -2*rz*cos(phi_z(k)), rz^2];
    a(k,:) = [1, -2*rz*cos(phi_p(k)), rp^2];
    Htmp = freqz( b(k,:), a(k,:), L );
    H(k,:) = Htmp(:).';
end

% add real-valued zero at DC
[H0,w] = freqz( [1, -rz], 1, L );
H = [ H0(:).'; H ];

Htot = prod(H);
semilogx( w/pi+1e-5, abs(Htot) );
axis([1e-5,.01,0,.0005])


Multiplying signals in discrete-time vs continuous-time Given two discrete-time signals $a[n]$, $b[n]$ and its product $c[n]=a[n] b[n]$. The ideally interpolated, continuous-time version of $c[n]$ is
\begin{align}
c_1(t)&=\sum_{n=-\infty}^{\infty} a[n] b[n] \operatorname{sinc}\left(\frac{t-n T}{T}\right) .
\end{align}
Now I would like to perform this multiplication in continuous-time using the interpolated versions of $a[n]$ and $b[n]$:
\begin{align}
c_2(t) &= \left(\sum_{n_1=-\infty}^{\infty} a[n_1] \operatorname{sinc}\left(\frac{t-n_1 T}{T}\right)\right) \left(\sum_{n_2=-\infty}^{\infty} b[n_2] \operatorname{sinc}\left(\frac{t-n_2 T}{T}\right)\right) \\
&= \sum_{n_1=-\infty}^{\infty}\sum_{n_2=-\infty}^{\infty} a[n_1] b[n_2] \operatorname{sinc}\left(\frac{t-n_1 T}{T}\right)   \operatorname{sinc}\left(\frac{t-n_2 T}{T}\right) .
\end{align}
Suddenly I end up with a double-sum and the cross terms $n_1\neq n_2$ are nonzero.
Why?
What do I need to do to arrive in a form similar as $c_1(t)$ with just one sum?

You reached a puzzling conclusion about $c_1(t) = c_2(t)$, and wonder whether you made a mistake in deriving them, or if the equality is indeed correct then how to explain it, perhaps by explicitly deriving one from the other.
I cannot tell whether it's possible to explicitly manipulate the double summation in $c_2(t)$ so as to convert it into the single summation of $c_1(t)$. It may be possible, I haven't tried. But sometimes it's not possible to explicitly show it, and instead, we have to rely on indirect evidence to confirm the equality and use the equality for our advantage if possible.
One such example is the Fourier transform of the constant $1$ :
$$ \mathcal{F}\{1\} = 2\pi \delta(\omega) \tag{0} $$
The equality in Eq.0 is not derived by explicitly evaluating the forward Fourier integral, but instead, deduced from  the duality property of CTFT, and given that Fourier transform of $\delta(t)$ is $1$.
At the end of the analysis, we conclude that the validity of the equality $c_1(t) = c_2(t)$ is a consequence of the Nyquist sampling theorem;i.e., the truth of the equality is imposed by the sampling theorem, rather than a result of explicit algebraic manipulations of $c_2(t)$ into $c_1(t)$ or vice versa. And indeed, this is a useful side application of the theorem to prove that some equation is true which is very hard, if not impossible, to do so otherwise.
Let me show you, therefore, an indirect way of imposing the equality.
Let all signals $a(t),b(t)$, and $c(t)=a(t)b(t)$ are sufficiently bandlimited so that we can avoid aliasing.
Observe impulse train modulation relation:
$$ x_s(t) = x(t) \cdot \delta_T(t) ~ \cdot \cdot \cdot ~ \delta_T(t) = x(t) \cdot \delta_T(t) \tag{1}$$
where  $ ~\delta_T(t) = \sum_n \delta(t-nT)$.
Also observe the interpolation relation :
$$ \left( x(t) \cdot \delta_T(t) \right) \star h(t) = x(t) \tag {2}$$
where $~h(t) = \text{sinc}(t/T) ~$ is the ideal lowpass brickwall interpolation filter.
We will use Eqs. 1 & 2 to derive alternative but equivalent expressions for the samples $c_s(t)$ of $c(t)$ and achieve the single and double sum versions $c_1(t)$ and $c_2(t)$.
The Double Sum :
$$
\begin{align}
c(t) &= c_s(t) \star h(t) \\
&= \left( c(t) \cdot \delta_T(t) \right) \star h(t) \\
&= \left( a(t) \cdot b(t) \cdot \delta_T(t) \right) \star h(t) &\text{Multiply then sample}\\
&= a(t)b(t)  &\text{by Eq.2} \\
&= \left( \sum_n a[n] h(t-nT) \right) \left( \sum_m b[m] h(t-mT) \right) \\
&= \sum_n \sum_m a[n] b[n] ~ h(t-nT) ~ h(t-mT) \\
\end{align}
$$
The Single Sum :
$$
\begin{align}
c(t) &= c_s(t) \star h(t) \\
&= \left( a(t) \cdot b(t) \cdot \delta_T(t) \right) \star h(t) \\
&= \left( a(t) \cdot \delta_T(t) \cdot b(t) \cdot \delta(t) \right) \star h(t) &\text{by Eq.1} \\
&= \left( a_s(t) \cdot b_s(t) \right) \star h(t) &\text{Sample then multiply} \\
&= \left( \sum_n a[n] \delta(t-nT) \right) \left( \sum_m b[m] \delta(t-mT) \right) \star h(t) \\
&= \left( \sum_n \sum_m a[n]b[m] \delta((m-n)T) \cdot \delta(t-mT) \right) \star h(t) \\
&= \sum_n a[n] \left( \sum_m b[m] \left[ \delta((m-n)T)\delta(t-mT) \star h(t) \right] \right) \\
&= \sum_n a[n] \left( \sum_m b[m] \delta((m-n)T)) h(t-mT) \right) \\
&= \sum_n a[n]b[n] h(t-nT)
\end{align}
$$
The sampling relations are  :
$$a_s(t) = a(t) \delta_T(t) = \sum_n a[n] \delta(t-nT) \tag{3}$$
$$b_s(t) = b(t) \delta_T(t) = \sum_m b[m] \delta(t-mT) \tag{4}$$
$$c_s(t) = c(t) \delta_T(t) = \sum_n c[n] \delta(t-nT) \tag{5}$$
We can get back $a(t),b(t)$,and $c(t)$ by ideal bandlimited interpolation of their samples $a[n],b[n]$,and $c[n]$ :
$$a(t) = a_s(t) \star h(t)= \sum_n a[n] \text{sinc}((t-nT)/T) \tag{6}$$
$$b(t) = b_s(t) \star h(t)= \sum_m b[m] \text{sinc}((t-mT)/T) \tag{7}$$
$$c(t) = c_s(t) \star h(t)= \sum_n c[n] \text{sinc}((t-nT)/T) \tag{8}$$


What does "Lag" mean in terms of Cross Correlation? I am working on a program using the Essentia library for audio analysis. I've implemented the Cross Correlation algorithm (here: http://essentia.upf.edu/documentation/reference/std_CrossCorrelation.html), but I'm not sure if I really understand the output. Sadly, barely any of the description on that page makes sense to me. I do not have a math background, and only a cursory understanding of frequency analysis / signal processing.

So... can anyone explain what the optional parameters (minLag, maxLag) do here?

using the simplest definition of the cross-correlation:

$$ R_{xy}[k] \triangleq \sum_n x[n]\,y[n+k] $$

the "lag" is the displacement $k$.  (i am being deliberately vague about the limits to the summation.)  in the correlation $x[n]$ is lagging behind $y[n]$ by $k$ sample periods.


Why would you use DCT as oppose to DFT for detecting oscillation in a signal? For my project, I needed to detect the primary oscillation frequency in a signal data. As I googled around, it turns out that DCT is the best tool to find this (http://www.mech.pku.edu.cn/robot/teacher/wangjiandong/WangHuangLu_ImprovedDCT_CEP_2013.pdf)

My question is: why would you base the analysis on DCT, and not DFT? What would happen if I were to apply the same analysis but with DFT? I can't accept that DCT is "faster" since for the practice of detecting oscillation, even in real time, both are "fast enough" and so speed and computation complexity should not be a factor into deciding which to pick!

Just keepin' it $\Re$eal.

The DFT records the timing offset of a particular sinusoid frequency as the phase of the complex coefficient.  Since the DFT of a real signal is conjugate-symmetric, the corresponding negative frequency coefficient gets the same magnitude but negative phase.  FFT routines frequently skip calculating the negative spectrum for real signals.
This leaves NFFT/2+1 non-redundant complex coefficients, however the DC and Nyquist bins are completely real.  It would be possible to combine the DC and Nyquist bin together, making zero growth in data size (for even length signals anyway).

But this kind of flipping between Complex and Real structure may not be convenient for analysis or computation.

The DCT captures essentially the same information about the real signal, except the timing of the sinusoid is encoded as the sum of cosine waves.
The DCT has the advantage of representing a unitary transform without mixing real and complex numbers.  The DFT is also unitary, but it maps $\mathbb{R^N} \rightarrow \mathbb{C^N}$.


Can the instantaneous frequency be always derived from an analytic signal? Instantaneous frequency can be defined as a derivative of an instantaneous phase of an analytic signal which can be nicely seen in practice in this example from Scipy's documentation. But it seems like it does not always work like this. I played with the code from the example and obtained varied results like for this sine wave with frequency gradually changing from 2 to 12Hz:

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import hilbert, chirp

duration = 1.0
fs = 400.0
samples = int(fs*duration)
t = np.arange(samples) / fs
w = 2*np.pi*(2 + 10*t)
signal = np.sin(w*t)
analytic_signal = hilbert(signal)
amplitude_envelope = np.abs(analytic_signal)
instantaneous_phase = np.unwrap(np.angle(analytic_signal))
instantaneous_frequency = (np.diff(instantaneous_phase) /
                           (2.0*np.pi) * fs)

This is the resulting figure with the original frequency plotted in grey on the bottom subplot:

fig = plt.figure()
ax0 = fig.add_subplot(211)
ax0.plot(t, signal, label='signal')
ax0.plot(t, amplitude_envelope, label='envelope')
ax0.set_xlabel("time in seconds")
ax0.legend()
ax1 = fig.add_subplot(212)
ax1.plot(t, 0.5*w/np.pi, lw=2, color='gray')
ax1.plot(t[1:], instantaneous_frequency)
ax1.set_xlabel("time in seconds")
fig.show()

In the bottom plot the blue line and the grey line go in the same direction, but there is a huge discrepancy between the two and I do not mean the ever-present edge effect. The calculated instantaneous frequency is rising at the rate that is almost a double of the actual reaching values above 20Hz at the end.

I tried it for other sine waves with constant amplitudes and varied frequencies getting similar results. Only when I used a signal with a constant frequency the blue line was matching the grey line on the bottom plot. My question here is: are there any conditions a signal must meet for the definition of the instantaneous frequency to hold true, so it can be derived from the instantaneous phase of the analytic signal?

This is because the phase term is changing with time at the rate of $t^2$ for the frequency ramp, and since the instantanous frequency is the derivative of the phase your result is as expected.

$m(t) = cos(2\pi f(t) t)$

$\phi = 2\pi f(t) t$

$f(t) = Kt$ for a linear frequency ramp

thus phase versus time changes at the rate of $t^2$

$\phi = 2\pi Kt^2$

and the derivative of the phase (which is the instantaneous radian frequency) is:

$\frac{d\phi}{dt} = 4\pi Kt$

(Divide this by $2\pi$ to get instaneous frequency, but the source of the 2x discrepancy should now be very clear).

See the answer to this DSP Challenge that I posted that further explains this:

Simulation of a Frequency ramp


What is the difference between MSE and PSNR of images? I have some doubts.

What is the difference between quality parameters MSE and PSNR?
Which is used when?
Is there any advantage of PSNR over MSE or vice versa?

Also I have heard people saying that MSE and PSNR should never be used to check image quality. Instead one should use UQI and SSIM. I read about UQI and SSIM and came to know they are close to human visual system but I read a technical report which says that SSIM is sometimes better than MSE and sometimes MSE is better. Here is the link
http://www.cs.uregina.ca/Research/Techreports/2008-02.pdf

Now I am confused which one to use(among MSE and SSIM) and why?

Also it is said that SSIM is upgraded version of UQI. Then what's the difference between these two?
Why do people still use UQI irrespective of the fact that SSIM is there?

These metrics are used to compute certain signal properties. They are roughly divided into 2 categories as 1-objective, 2-subjective. Objective criteria are too mechanical and do not take any human perception weighting. Subjective criteria are either model based computations or directly opinion based measurements. Which is to use when is a matter of application.


Similarity in signal plots but dissimilarity in sensed content I have two audio signals that look incredibly similar when plotted. However, when I play them, they sound quite different. The only difference I can even see in their plots is that one signal looks darker than the other. Why would this be happening?

FYI: The signal on the top is one that has been reconstructed from a reconstruction algorithm and the one in the bottom is the original one.

Here's some further data about how the signal has been produced. The original signal is a sample audio data I downloaded off the internet. The purpose of my experiment is to produce a close enough signal from only the amplitude spectrum and without the phase information. For similar, but slower algorithms, one can look at the Griffin and Lim algorithm and also the single pass spectrogram inversion paper.

In my case, I'm getting a waveform very close to the original signal, but there's a significant "noise" when listening to the output waveform.

So in short, the question is: waveforms looking so close to each other sound very different which is not very intuitive. Also, when signals look so similar, it is difficult to discern that there's any kind of noise involved. So what kind of hidden noise is this?

One last piece of information: the correlation function gives 0.69 as the output when comparing these two signals. Not sure if that's useful

First, the original $x$ and reconstructed  $\hat{x}$ signals have a peak amplitude around $0.35$. Their difference peaks above $0.40$. That happens between two time-shifted like-alike signals.

Second, the difference seems to relate to the amplitude. Processing maybe be non-linear, and it could be interesting to look at relative differences like $2(x -\hat{x})/(x +\hat{x})$ as well.

To enhance understanding, what would  3D plots $(x[n],x[m],x[n]-x[m])$ and $(x[n],x[m],2\frac{x[n]-x[m]}{x[n]+x[m]})$ like like?


Approximation using a Fourier transform with low pass filter I need to approximate a function f, but I cannot do so with frequencies that exceed 1kHz

What is the best approximation I can get? Is taking the Fourier transform then zeroing any term above 1kHz the best approximation? Or can I fiddle with the lower order terms and get a better fit?

This is a made up scenario, but I have to prove the same concept with Walsh transforms. I am fairly certain that the lower order terms form the best approximation from random twiddling and hill climbing searches, but I need proof.

I believe the proof is something very similar to a least squares regression proof, but I can't get it. Has this problem been solved before? At least in the Fourier domain?

Depends on what you mean by 'best'.  In terms of least squares, then yes, just knocking off the terms above 1kHz will give the closest approximation.  I guess one way to prove it would be:
$$
     x(t) = \int_{-\infty}^{\infty} X(f)\exp\left( j2\pi ft \right) dt\\
\hat x(t) = \int_{-1k}^{1k} X(f)        \exp\left( j2\pi ft \right) dt
$$
leaving a residual error of
$$
e(t) = \int_{|f| > 1k} X(f)\exp\left( j2\pi ft \right) dt .
$$

Compare this to a 'fiddled' version, using $\tilde X(f)$ for the coefficients:
$$
\hat x(t) = \int_{-1k}^{1k} \tilde X(f) \exp\left( j2\pi ft \right) dt,
$$
which has a residual error of
$$
e(t) = \int_{|f| \le 1k} \left(X(f) - \tilde X(f)\right) \exp\left( j2\pi ft \right) dt +
\int_{|f| > 1k} X(f)\exp\left( j2\pi ft \right) dt ,
$$
which is greater than the error in the first instance unless $\tilde X(f) = X(f) \; \forall f$.

Something like that ought to do it, anyway.

EDIT: It is possible to prove that a projection onto any orthogonal set of functions (doesn't need to be the Fourier set) gives the least squares fit:

$$
\int f_k(t) f_m(t) dt = \cases{
  0, & $k \ne m$ \\
  1/\lambda_k, & $k = m$}
$$
for any orthogonal basis set $\{f_k\}$ (I'm not bothering with weighting functions), where
$$
\lambda_k = \frac{1}{\int f^2(t) dt}
$$
The projection of $g(t)$ onto the set is given by
$$
c_k = \lambda_k \int g(t)f_k(t) dt.
$$
The least squares fit of a set of orthogonal basis function is the minimisation of
$$
\int \left[g(t) - \sum_k d_k f_k(t)\right]^2 dt.
$$
Minimise by setting the derivative (wrt $d_m$) to zero:
$$
2\int\left(g(t) - \sum_kd_kf_k(t)\right)\left(-f_m(t)\right) dt = 0 \\
-\int f_m(t)g(t)dt + \int f_m(t)\sum_kd_kf_k(t)dt = 0 \\
-\int f_m(t)g(t)dt + \sum_kd_k \int f_m(t)f_k(t)dt = 0 \\
-\int f_m(t)g(t)dt + \sum_kd_k \frac{1}{\lambda_m} = 0
$$
which rearranges to give the same value as for $c_k$, thus the least squares fit is exactly the straightforward projection onto the orthogonal set.


Retaining filter state for processing data in chunks I am have data in 1 minute chunks but I would like perform filtering on these as if the data was one contiguous signal. How might I hold the filters state between processing data chunks to achieve this.

I am using an FIR in MATLAB with the filter() function.

You can use the filter function of matlab with initial conditions with the following syntax: (excerpted from Matlab 6.0 command window help)

[Y,Zf] = FILTER(B,A,X,Zi) gives access to initial and final conditions, Zi
and Zf, of the delays.  Zi is a vector of length MAX(LENGTH(A),LENGTH(B))-1
or an array of such vectors, one for each column of X.

Hope this helps.


Impulse response of an infinite-length system I have a system which is expressed as following:
$$y[n] = x[n] - \alpha y[n-N]$$ where $y[n]$ is the output and $x[n]$ is the input.
When I take the z-transform of both sides, I've found the transfer function as $$H(z)=\frac{1}{1+\alpha z^{-N}}$$.
However, I need to find the impulse response of this system but I couldn't find $h[n]$. Can anybody help me how to solve this problem?

Cconsider the following pair of signals and their transforms:
If $g[n]$ has a Z-transform $G(z)$ then the expanded (zero stuffed) signal $h[n]$ has the Z-transform $H(z) = G(z^N)$, where
$$ h[n] = \begin{cases} {g[\frac{n}{N}] ~~~,~~~ n = m N  \\ 0 ~~~~~~~~~,~~~ \text{otherwise}  \tag{1}}\end{cases} $$
Looking at the given Z-transform :
$$H(z) = \frac{ 1}{ 1 - \alpha z^{-N} }$$
it can be seen that if $H(z) = G(z^N)$ where
$$G(z)= \frac{ 1}{ 1 - \alpha z^{-1} }$$
and $g[n] = \alpha^n u[n] $.
Then you can simply conlcude that :
$$ h[n] = \alpha^{n/N} ~u[n/N] =  \begin{cases} {\alpha^{n/N}  ~~ ,~~ n = m N, m = 0,1,..  \\ 0 ~~~~~~~~~,~~ \text{otherwise} }\end{cases} $$


Disadvantage of high downsampling rates with CIC filters For a lock-in amplifier I am downsampling the demodulated baseband signal  prior to the final low pass filtering. This lock-in amplifier is implemented on an FPGA (myRIO of National Instruments).

Therefore I use CIC filters. I've reading a lot in Hogenauers paper and in "Understanding digial signal processing", by Richard Lyons. My downsampling rates are between 512 and 2048 (always powers of 2) as I want to decimate my input signal down to less than 100 Hz (e.g. 48kHz sampling will be decimated to 93Hz).

The gain of the CIC output is scaled to the CICs input by a simple bitshift afterwards.

Is there any disadvantage using even higher decimation rates?

CIC filters don't have very steep roll-off, or a very flat passband. That's the downside of using this very special class of filters.

Other than that, you should be fine. However, as with every filter design, you will have to design, and verify, yourself to ensure the filter matches your requirements.

However, CICs are very often used as runtime-adjustable filters; and that works pretty well if you have e.g. high-quality Nyquist-M (especially halfband) FIRs. An example of this is the Ettus USRP signal processing DDC:

   // CIC decimator  24 bit I/O
   cic_strober cic_strober(.clock(clk),.reset(rst),.enable(ddc_enb),.rate(cic_decim_rate),
               .strobe_fast(1),.strobe_slow(strobe_cic) );

   cic_decim #(.bw(WIDTH))
     decim_i (.clock(clk),.reset(rst),.enable(ddc_enb),
          .rate(cic_decim_rate),.strobe_in(1'b1),.strobe_out(strobe_cic),
          .signal_in(i_cordic_clip),.signal_out(i_cic));

   cic_decim #(.bw(WIDTH))
     decim_q (.clock(clk),.reset(rst),.enable(ddc_enb),
          .rate(cic_decim_rate),.strobe_in(1'b1),.strobe_out(strobe_cic),
          .signal_in(q_cordic_clip),.signal_out(q_cic));

   // First (small) halfband  24 bit I/O
   small_hb_dec #(.WIDTH(WIDTH)) small_hb_i
     (.clk(clk),.rst(rst),.bypass(~enable_hb1),.run(ddc_enb),
      .stb_in(strobe_cic),.data_in(i_cic),.stb_out(strobe_hb1),.data_out(i_hb1));

   small_hb_dec #(.WIDTH(WIDTH)) small_hb_q
     (.clk(clk),.rst(rst),.bypass(~enable_hb1),.run(ddc_enb),
      .stb_in(strobe_cic),.data_in(q_cic),.stb_out(),.data_out(q_hb1));

   // Second (large) halfband  24 bit I/O
   wire [8:0]  cpi_hb = enable_hb1 ? {cic_decim_rate,1'b0} : {1'b0,cic_decim_rate};
   hb_dec #(.WIDTH(WIDTH)) hb_i
     (.clk(clk),.rst(rst),.bypass(~enable_hb2),.run(ddc_enb),.cpi(cpi_hb),
      .stb_in(strobe_hb1),.data_in(i_hb1),.stb_out(strobe_hb2),.data_out(i_hb2));

   hb_dec #(.WIDTH(WIDTH)) hb_q
     (.clk(clk),.rst(rst),.bypass(~enable_hb2),.run(ddc_enb),.cpi(cpi_hb),
      .stb_in(strobe_hb1),.data_in(q_hb1),.stb_out(),.data_out(q_hb2));

The job of this chain is to decimate a 100MS/s signal down, adjustably, down to rates of 1/512.

Generally, if you just want to use powers of two, you can just use the same halfband filter multiple times; if you're doing it elegantly, you can probably even multiplex the same hardware for multiple iterations.

By the way, coming from a background where sample rates are in the MS/s, I wonder why you need such a low-rate decimation in hardware – every CPU can do that, nowadays, without breaking a sweat.


Help me understanding this source code of a Phaser I tried to understand the following source code of a Phaser:

#define SR (44100.f)  //sample rate
#define F_PI (3.14159f)

 class Phaser{
 public:
     Phaser()  //initialise to some usefull defaults...
      : _fb( .7f )
      , _lfoPhase( 0.f )
      , _depth( 1.f )
      , _zm1( 0.f )
  {
    Range( 440.f, 1600.f );
    Rate( .5f );
  }

  void Range( float fMin, float fMax ){ // Hz
      _dmin = fMin / (SR/2.f);
      _dmax = fMax / (SR/2.f);
  }

  void Rate( float rate ){ // cps
      _lfoInc = 2.f * F_PI * (rate / SR);
  }

  void Feedback( float fb ){ // 0 -> <1.
      _fb = fb;
  }

  void Depth( float depth ){  // 0 -> 1.
      _depth = depth;
  }

  float Update( float inSamp ){
      //calculate and update phaser sweep lfo...
      float d  = _dmin + (_dmax-_dmin) * ((sin( _lfoPhase ) +  1.f)/2.f);
      _lfoPhase += _lfoInc;
      if( _lfoPhase >= F_PI * 2.f )
          _lfoPhase -= F_PI * 2.f;

      //update filter coeffs
      for( int i=0; i<6; i++ )
          _alps[i].Delay( d );

      //calculate output
      float y =     _alps[0].Update(
                   _alps[1].Update(
                    _alps[2].Update(
                     _alps[3].Update(
                      _alps[4].Update(
                       _alps[5].Update( inSamp + _zm1 * _fb ))))));
      _zm1 = y;

      return inSamp + y * _depth;
  }
 private:
     class AllpassDelay{
     public:
      AllpassDelay()
          : _a1( 0.f )
          , _zm1( 0.f )
          {}

      void Delay( float delay ){ //sample delay time
          _a1 = (1.f - delay) / (1.f + delay);
      }

      float Update( float inSamp ){
          float y = inSamp * -_a1 + _zm1;
          _zm1 = y * _a1 + inSamp;

          return y;
      }
  private:
      float _a1, _zm1;
  };

  AllpassDelay _alps[6];

  float _dmin, _dmax; //range
  float _fb; //feedback
  float _lfoPhase;
  float _lfoInc;
  float _depth;

  float _zm1;

};

I understood most of it but i can´t find out what the purpose of the _zm1 variables is. There are two of them one in the AllpassDelay class, which are never used and one in the Phaser class which modifies the feedback.
Can someone explain me what this does?
Thanks

_zm1 is the internal state variable of the allpass. The allpass is an IIR filter and hence it's current output is dependent on past inputs and outputs. This is the way to keep track of this.

In this particular case it's implemented using a set of coupled difference equation in so-called "Transposed Form II" filter topology. So $zm_1[n] = y[n]*a_1+x[n]*b_1 $ . For a first order allpass we have $b_1 = 1$


Improving template with instantaneous amplitude Given some recording (containing both the signal of interest and some Gaussian noise), can we improve our template (or model) of the signal using the analytic signal of the recording?

For example, a signal (red) with Gaussian noise results in some noisy recording (black). The instantaneous amplitude of the recording (black) deviates away from the true signal's instantaneous amplitude (red) because of this noise:

This got me thinking: Since we know the instantaneous phase of the true signal will be linear (when unwrapped), we can reconstruct this from the (noisier) instantaneous phase of the recording, say using OLS. For argument's sake, let's assume a fixed unit amplitude (even though in this example the amplitude is actually modulated at 2 Hz) and reconstruct the recordings either from the raw instantaneous phase (i.e., with noise) or from the reconstructed instantaneous phase (i.e., hopefully without noise) -- and compare these to the reconstructed signal from the instantaneous phase of the true signal (i.e., without noise).

Clearly, the reconstructed (OLS) instantaneous phase better approximates the true signal's instantaneous phase than the instantaneous phase we yield from the raw recording (error = true signal – reconstruction; figure  below).

Finally, to see how close the reconstructions (either using raw or OLS instantaneous phase) approximate the true signal, I've computed these using the raw instantaneous amplitude for both. While (again) the OLS does better (below), there is still so much noise (because the instantaneous amplitude is computed from the raw recording).

Is there a way to better approximate the true signal (e.g., by low-pass filtering the instantaneous amplitude)?

Since I don't really know what you're using software-wise, I just whipped up a really short demo in GNU Radio Companion. This uses a PLL in a live processing demo:

If you're inclined to play around with it, here is the flow graph ready to be run if you've got a running installation of GNU Radio. The flow graph looks like this:

The left part up to the conversion to complex (which actually is a bit superfluous) is just the generation of a noise cosine.

The "PLL Freq Det" Block does pretty much what you describe: it compares the phase of the signal it observes with an internally generated oscillator, and generates a "phase signal" that is used to tell the internal frequency reference "how fast" it should oscillate. That signal is just a phase difference to be compensated, and needs to be internally filtered (the loop filter does that) before being applied to the internal oscillator.

At the same time, that signal is actually also a phase per sample - i.e., a frequency. So, you get a frequency estimate out of this.

Now, this is really practical for things like radio receivers that need to update their estimate of the received signal's frequency while operating, but in your case, you can do all this much more accurate, since you can actually work on the full signal.

A very popular (albeit, seldom optimal) estimator for the spectrum of a signal (and you're looking for which frequencies are making up your signal, so that's a question about the spectrum of your signal) is in fact the (magnitude of the) discrete fourier transform (DFT), often implemented by an FFT.

If you, for example, have your data in a Python program as vector signal, the following (untested) code should do interesting things to your signal:

import numpy
from matplotlib import pyplot

## Replace with your own samples and rate
samples = …
sampling_rate = 32e3

## Do a fourier transform of the whole data set
fourier_transform = numpy.fft.fft(samples)
positive_freq_components = fourier_transform[0:len(fourier_transform) / 2]
positive_freq_components_mag_squared = numpy.abs(positive_freq_components)**2
frequencies = numpy.linspace(0, sampling_rate / 2, len(fourier_transform) / 2)

## Plot that stuff!
pyplot.plot(frequencies, positive_freq_components_mag_squared)
pyplot.xlabel("Frequency (Hz)")
pyplot.ylabel("Digital Power")
pyplot.tight_layout()
pyplot.show()

This should give you a plot of power over frequency. Your noise will probably be relatively white, i.e. it will not have a lot of concentration in frequency domain. The sine, on the other hand, is very periodic and thus has all its energy at one frequency.

Note that I said "seldom optimal" about the DFT as spectrum estimator: It gives you the same frequency resolution all over your bandwidth, and you don't care about most frequencies. The (usable) frequency resolution depends linearly on the observation length, and you can't just get arbitrary frequencies by looking for peaks – all frequencies represented are sampling rate / length of DFT.

Other, "cooler", methods can give you a spectrum estimate that you can evaluate at arbitrary frequencies (such as MUSIC¹), or directly an estimate for the "strongest" frequency (or the strongest N frequencies) (such as ESPRIT¹). However, using these often requires a few strong assumptions, so I won't go deeply into that: If you want to know more about parametric estimators for fundamental frequencies, feel encouraged to research and ask a new question – there's folks on here way more experienced with such estimation tasks than I am.


How to find inverse z transform Suppose
$$Y(z) = \frac{\frac 12 z + 1}{z+\frac 12} \cdot \frac{z}{z-\frac12}\text.$$
According to Wolfram Alpha the inverse transform is, $2^{-n - 2} \cdot(5 - 3 \cdot (-1)^n)$. However, I cannot show why this is true. Could someone please guide me?

Using binomial formulas one easily gets
$$
Y(z)=\frac{\frac12+z^{-1}}{1-\frac14z^{-2}}
=(\tfrac12+z^{-1})\sum_{k=0}^\infty 2^{-2k}z^{-2k}.
$$
Now compare to the coefficients in $Y(z)=\sum_{n=0}^\infty y_nz^{-n}$ for $n=2k$ and $n=2k+1$, even and odd indices.


How to find H(z) and H(k) from a given causal function Consider the causal function,
$y[k] = 2x[k] - 40x[k - 1] + 10y[k - 1]$  $16y[k - 2]$;
where $y[k]$ is the output and $x[k]$ is the input. Assume that the system is initially at rest.
Please someone help me find $H(z)$ and $H(k)$ from given equation. I have an assignment due to tomorrow

It is a Homework problem, so only Hint will be provided. You need to show some work of your own following the hints.

First of all, take $\mathcal Z$-transform of the equation you have. Take care of the delays in time-domain and apply appropriate property of $\mathcal Z$-transform.

Then, try to get the transfer function $H(z)$ from there.

Once you have $H(z)$, figure out a way to convert $H(z)$ into $H(e^{j\omega})$ and then try to proceed for $H[k]$.

Show your work by making updates in your question and we will be clarify the doubts and be able to help you understand.