What is the physical meaning of coherence bandwidth in mmWave massive MIMO channel? I am doing some research about millimeter wave massive MIMO channel now.
We first take account about static MIMO channel, cause we only have SISO measurement system now.
We use scanner to move our antenna to simulate it as array antenna to measure the static channel.
Afterwards, using signal processing method to simulate static MIMO channel.
There are some paper that doing the similar research, they analysis the coherence bandwidth at different antenna element. But I can't realize what the purpose of studying that issue. Is that related to some other MIMO  techniques?

As far as I know, coherence bandwidth is corresponding to multi-path delay spread, and is an indicator to distinguish whether the channel is flat-fading or not.

In the other point side of view, I also study multi-path MIMO channel model. The model is represent as below

It shows that multi-path would cause difference phase shift owing to the steering array. However, I can't realize what the model is associated with coherence bandwidth at each element.

Can someone please help me or provide some idea. Thanks in advance

reference ： R. W. Heath, N. González-Prelcic, S. Rangan, W. Roh and A. M. Sayeed, "An Overview of Signal Processing Techniques for Millimeter Wave MIMO Systems," in IEEE Journal of Selected Topics in Signal Processing, vol. 10, no. 3, pp. 436-453, April 2016.

You are right, coherence bandwidth is the frequency domain counterpart of delay spread.

However, to achieve diversity across antennas the spacing between the antennas is important "relative to the environment".Let me explain that a bit more. In a user mobile phone you would typically find the order of wavelengths seperation between antennas sufficient to get independent fading accross them. However, the same is not true at the base station. Typically at the base station the seperation between the antennas has to be a 10s of order of wavelength. Why so? Because a user equipment near to ground sees a much richer scattering environment, which are enough to cause phase changes within a wavelength. However at the base station (usually placed high up in the air) the scattering environmental is not that rich and hence more separation is required.

So, to study the effect of independent fading with antenna separations one can use the methodology to see if other antennas are seeing the same coherence bandwidths or not. It is all part of channel profiling whenever new frequencies are deployed.


How to assess perceptual volume from FFT So far I could estimate signal volume from FFT using amplitude information, then adding all amplitudes together and divide this number by the number of frequencies. But I think this method does not apply for estimation of subjective perceptual volume.

If you equate perceptual volume to loudness, there are international standards.

if you look at
https://www.mathworks.com/matlabcentral/fileexchange/46819-a-weighting-filter-with-matlab-implementation

This is an A-weighted Matlab FFT based implementation.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                  A-weighting Filter                  %
%              with MATLAB Implementation              %
%                                                      %
% Author: M.Sc. Eng. Hristo Zhivomirov        06/01/14 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function xA = filterA(x, fs)

% function: xA = filterA(x, fs)
% x - signal in the time domain
% fs - sampling frequency, Hz
% xA - filtered signal in the time domain

% determine the signal size
sz = size(x);

% represent x as column-vector
x = x(:);

% signal length
xlen = length(x);

% number of unique points
NumUniquePts = ceil((xlen+1)/2);

% FFT
X = fft(x);

% fft is symmetric, throw away second half
X = X(1:NumUniquePts);

% frequency vector with NumUniquePts points
f = (0:NumUniquePts-1)*fs/xlen;

% A-weighting filter coefficients
c1 = 3.5041384e16;
c2 = 20.598997^2;
c3 = 107.65265^2;
c4 = 737.86223^2;
c5 = 12194.217^2;

% evaluate the A-weighting filter in the frequency domain
f = f.^2;
num = c1*f.^4;
den = ((c2+f).^2) .* (c3+f) .* (c4+f) .* ((c5+f).^2);
A = num./den;
A = A(:);

% filtering in the frequency domain
XA = X.*A;

% reconstruct the whole spectrum
if rem(xlen, 2)                     % odd xlen excludes the Nyquist point
    XA = [XA; conj(XA(end:-1:2))];
else                                % even xlen includes the Nyquist point
    XA = [XA; conj(XA(end-1:-1:2))];
end

% IFFT
xA = real(ifft(XA));

% represent the filtered signal in the form of the original one
xA = reshape(xA, sz);

end

The Mathwork will also sell you their Audio toolbox that has many types of weightings.

I find these standards a bit confusing because you don't know if they correct for the effective noise bandwidth of proportional bandwidth filters and if an FFT implementation should take that into account.

I don't vouch for the Matlab code above.  It seems reasonable.

The Wikipedia article:

https://en.wikipedia.org/wiki/A-weighting

also gives the analog transfer functions, which can be converted to digital filters, but actual digital filters seem harder to find.

If your ultimate purpose is to use a standard where a professional signature is required, or one is exposed to liability, I would be inclined to use the Mathwork's product.


Efficient computation of dot product of convolutions I have two FIR filters $f_{1},f_{2}$ and two short time-slices of signals $s_{1},s_{2}$. I need to compute: $$\left<s_{1}*f_{1},\ s_{2}*f_{2}\right>$$

Which I compute using six FFTs: each $s_{i}*f_{i}$ is computed by pointwise multiplication of the (padded) FFT vector of $s_{i}$ by the corresponding FFT of $f_{i}$, then an IFFT is computed and the result is truncated (this truncation is the reason computing the dot product in the frequency domain does not work).

The $f_{i}$ are known ahead of time, so this takes four FFTs at each iteration of my code.

Is there a way to get away with computing less? If the conjugate transpose of $f_{1}$ could be composed on $f_{2}$ and a FIR filter was the result, it would be posible to manage with just 2 FFTs with something like: $$IFFT\left(FFT(s_{2})\cdot FFT(f_{2}^T\circ f_{1})\right)$$

(where $FFT(f_{2}^T\circ f_{1})$ is computed ahead of time,) but this doesn't quite seem to work out.

A team member of mine asked a similar question (on another stackexchange community) but did not get helpful replies.

edit (details about truncation in my code):

To compute $\left<s*f\right>$ (supposing that $s$ is of length 20000 and $f$ of length 4000, though these numbers are made up)  both vectors are zero-padded to length $20000+4000-1$, FFTs are multiplied pointwise, IFFT is computed on the product, and the samples in indices $\left[3999, 19998\right]$ are taken from the result (this represents a closed range of indices on a vector indexed from 0.)

This gives us all of the "meaningful" result computable from the slice of the signal: all dot products between $f$ and a shift of the signal-slice $s$ such that $f$ is shifted onto a completely-known part of $s$ (this is done over successive time slices of a signal as it is recorded.) It's important to notice the result is of length 16000 (the length of $s$ minus the length of $f$).

Convolution is a specific type of linear operator. So in all generality, you can write your expression as $\langle F_1 s_1 | F_2 s_2 \rangle$ with two linear operators $F_{1,2}$. Using the rules of operators in bilinear forms, we can rewrite this as $\langle s_1|F_1^\dagger F_2 |s_2\rangle$ and that again as $\langle s_1|\left(F_1^\dagger F_2\right) s_2\rangle$.

In your case, the two operators have a pleasant property. They are toeplitz, and if we take your padding into account, even circulant. The product of two circulant operators is again circulant, as is the conjugate transpose. That means you can write the result of the operator product as a single convolution. Your convolution kernel is the correlation product $\mathrm{corr}({f_1}^*,f_2)$.

If you implement this using fast convolution, you need to pad this correlation product and the size of your kernel will increase to the sum of the sizes of $f_{1,2}$. Accordingly, the required padding for $s_2$ will also increase.


Why does Fourier Series have $\sin$ and $\cos$ Components While we look at Fourier Series there are both $\sin$ and $\cos$ components.But I think $\sin$ component is ony needed to describe wave.why there is also an $\cos$ component in Fourier Series?

$S_n(x)$=$A_0 + \sum_{n = 1}^{\infty} (A_n\cos(nx) + B_n\sin(nx))$

Note that the value of the series

$$\sum_{n=1}^{\infty}B_n\sin(nx)\tag{1}$$

for $x=0$ is always zero, regardless of the coefficients $B_n$ (assuming convergence). Furthermore, since the functions $\sin(nx)$ are odd, (1) can only represent odd functions. So you cannot represent general functions with the series (1). This is why you also need the cosine terms. As an alternative, you could as well use a representation with phases in the argument of the sine:

$$f(x)=A_0+\sum_{n=1}^{\infty}C_n\sin(nx+\phi_n)\tag{2}$$

In this case you don't need any cosine terms because of the phases $\phi_n$. The representation in (2) is equivalent to the series in your question using sine and cosine terms. You can see this by using the identity $\sin(\alpha+\beta)=\cos\alpha\sin\beta+\sin\alpha\cos\beta$:

$$\sum_{n=1}^{\infty}C_n\sin(nx+\phi_n)=\sum_{n=1}^{\infty}\left[C_n\sin\phi_n\cos(nx)+C_n\cos\phi_n\sin(nx)\right]$$

So the coefficients $A_n$ and $B_n$ in your question are

$$A_n=C_n\sin\phi_n\\
B_n=C_n\cos\phi_n$$

And the other way around you have

$$C_n=\sqrt{A_n^2+B_n^2}\\
\phi_n=\arctan\left(\frac{A_n}{B_n}\right)\quad(\pm\pi)$$


Do you recognize this filter? I'm trying to decipher the uncommented source code of SoundBox which is a little music tracker written in javascript.

There is an included synthesizer, which applies a filter on the output samples. I'm struggling to understand how it works.
You can find where the filter is applied in the source code here (line 251). However it is not commented.

Here is my interpretation of its definition.
First, the filter parameters:

a = cutoff/samplerate
f = 1.5 * sin(PI * a)
q = 1.0 - resonance

cutoff is in Hz. The samplerate is fixed at 44100.

then the filter process

low[i] += f * band[i-1]
high[i] = q * (x[i] - band[i-1]) - low[i]
band[i] += f * high[i]

y[i] = low[i]

where x is the input sample array, y the output and low/high/band some other variable whose purpose I don't fully understand yet.

This looks like some kind of IIR filter to me. Do you recognize it? Can you break it down for me? Explain it to me?

This is the low-pass filter version. One is supposed to be able to get a high-pass filter (or a band-pass) by changing the output to y[i] = high[i] (or band[i]).

I've made some experiments and it works pretty well, however the actual cutoff is a bit shifted and is not exactly equal to cutoff. I've been able to correct the problem by changing the definition of f to f = sin(2.0 * PI * a). Now the cutoff is exactly where it is supposed to be, but I don't know why. If you want to see nice pictures of that, you can see the issue I've opened on GitHub, but no response yet.

Thanks.

This is a digital state variable filter. Note that your pseudo-code of the filtering process is not correct. You replaced all += in the original code by =, which results in a filter that wouldn't work as desired.

The actual recursive equations implemented by the code are

$$\begin{align}l[n]&=l[n-1]+fb[n-1]\\
h[n]&=q\big(x[n]-b[n-1]\big)-l[n]\\
b[n]&=fh[n]+b[n-1]\end{align}\tag{1}$$

where $l[n]$, $h[n]$, and $b[n]$ are the low-pass, high-pass, and band-pass outputs, respectively. The equations $(1)$ implement a second-order IIR filter (actually three).


Sound absorption coefficients for frequencies higher than 8kHz and lower than 80Hz? I am building a matlab simulation that applies room compensation to a 3D audio rendering system, and I'm searching for the frequency dependent sound absorption coefficients of various materials in order to extract the reflection coefficients.

However all the tables I found provide absorption coefficients for different octave frequency bands that have center frequency from 125Hz to 4kHz or at most 8kHz, but in my case I need to consider frequencies up to 20kHz.

My question is why are tables for sound absorption coefficients limited to such frequency range and how should I consider the coefficients for higher frequencies(i.e. higher than 10kHz) and also for lower ones(i.e. lower than 80 Hz)?

I found the tables in the annex section of the book "Auralization" by Dr. Michael Vorländer, which uses the ISO 354:2003 standard in order to perform the measurements.

The only resource where I found tables for frequency values higher than 8kHz  was this paper: https://www.degruyter.com/downloadpdf/j/aoa.2013.38.issue-2/aoa-2013-0020/aoa-2013-0020.pdf where it is proposed an alternative method w.r.t. the previously mentioned standard in order to take the measurements from 50Hz to 50kHz, however these are provided just for very few materials.

The main reason why the tables don't have it, is that it's hard to measure. The measurement technique in ISO 354:2003 relies on measuring the difference in reverberation times in a reverberation rooms with and with/out a material sample.  At higher frequencies, the reverb time is dominated by air absorption and and the sound field becomes less and less diffuse which violates the basic assumption behind the measurement. It's not uncommon to measure absorption coefficients larger than 1. That's not good !

I wouldn't mingle air absorption and wall absorption. They are different effects and should be modelled separately. They are basically additive. At 20 kHz air absorption is around 50-ish dB/100m. The early reflections in a room have a travel distance of only a few meters so there is some attenuation, but there is still plenty of energy left.

Keep in mind that wall absorption measurement is NOT done at 8 kHz but it's averaged over the 8 kHz octave, so the data covers frequencies up to 11.5 kHz. If you need to go hire, you can often simply extrapolate to 16 kHz by looking at the difference between the 4 kHz and 8 kHz value.

Easiest is to extrapolate in dB for the reflection coefficient. Example: 4 kHz = 0.6, 8 kHz = 0.8. Reflection coefficients are 0.4 and 0.2 or -4dB and -7dB respectively. So it's 3 dB drop of reflected energy per octave. Extrapolation yields -10 dB for 16 kHz, which is a reflection coefficient of 0.1 or an absorption of 0.9.

Something like this
$$\alpha_{16} = 1-10^{0.1 \cdot (2 \cdot 10 \cdot \log_{10}(1-\alpha_8)-10 \cdot \log_{10}(1-\alpha_4))}$$


What is a good way to locate butterfly in an image I have some images of butterflies. I've downloaded them from here. A sample is below:

In each image, the most prominent object is the butterfly (with flowers, leaf etc in background). I usually go for contour/edge detection to locate features in image, but in this case these algorithms are not working as I expected.

I want to locate butterfly and bound it with a box, or better, compute its contour. Is there a way to do this automatically using python-opencv or some other tool?

The butterflies are really different and the backgrounds as well, then it makes it really difficult. I would try one of these ideas:

Train the Viola & Jones algorithm on your butterflies. It's a cascade of weak classifiers usually used for face detection (algorithm implemented in all the cameras), but it can be trained for other purposed.
The background seems to be pretty different of the butterflies. I would try to apply the thresholding by histereis according to the values of the pixels on the borders. Like that you select all the similar pixels that touch the border, and it will reduced a lot the region of interest.


What's the difference between a DSP and microcontroller? What are the core attributes that make a DSP and microcontroller different from each other?

I have in mind a machine vision device that needs to run image processing algorithms but also control some networking and hardware I/O. How good would the DSP be at controlling hardware output pins, and how well can microcontrollers run image processing algorithms?

a Digital Signal Processor is one that has, in its instruction set, some instructions and addressing modes that are optimized for processing digital signals.

usually these optimizations can be shown around what is needed to perform the dot-product needed for an FIR filter.

$$ y[n] = \sum\limits_{i=0}^{L-1} h[i]\,x[n-i] $$

to do this in, say, $L$ instructions, a DSP must be able to do in one instruction:

multiply $h[i]$ and $x[n-i]$ together.
accumulate that product into an existing sum.
fetch the next $h[i+1]$ and $x[n-i-1]$ in anticipation of the next multiply-accumulate.  since these are two numbers to fetch, a DSP will use something called a Harvard architecture that has at least two separate memory spaces for $h[i]$ and $x[n]$ so the DSP can fetch these two numbers simultaneously.
addressing $x[n]$ must be in a circular queue.  a DSP will perform the modulo (or "wrap around") arithmetic on the index or address of $x[n]$ necessary without additional instructions.
the result $y[n]$ will eventually go to an output DAC or fixed-point stream and there is some way to saturate the value of $y[n]$ against some $\pm$ maximum without additional instructions.  if the DSP is a fixed-point DSP, then this accumulator register will have width in bits that is the sum of the bitwidth for $h[i]$ and the bitwidth for $x[n]$ plus a few more bits on the left as "guard bits".

a general purpose CPU can do all of these, but not likely in a single instruction cycle and things like modulo arithmetic and saturation will need their own specific instructions in a general-purpose CPU.

a DSP may also have some instructions and an addressing mode that facilitates the operations of the Fast Fourier Transform (FFT).  this may include instructions necessary to perform an FFT "butterfly" in as few as four instruction cycles (or two instruction cycles if SIMD is operational).


Wavelets in time series predictions Wavelet transforms are linear in nature. Then how good they are for time series predictions. Can they analyze nonlinear relations in data well??

Wavelets can be useful feature detectors, and one of their standard advantage is to remove some trend or background effects, due to the vanishing moment property.

So first, this can act, sometimes, as a very low level non-linear effect remover.

Second, one of the successes of wavelets is the nonlinear approximation, where coefficients are selected upon their (highest) amplitude. The shrinkage/thresholding approach is nonlinear, as the shrinkage of a sum of two signals is not the same as the sum  of the shrinkage of two signals.

Third, more generally, wavelet coefficients are often used as the input of nonlinear predictors, for instance neural networks.


How to set threshold value for Hough Transform I am new to Computer Vision, I am writing my own code for Hough Transform I have found the following algorithm:

Algorithm:

• The cell (i,j) corresponds to the square associated
with parameter values (θj, ρi).

• Initialize all cells with value 0.

• For each foreground point (xk,yk) in the thresholded
edge image
– Let θj equal all the possible θ-values

• Solve for ρ using ρ=x cos θj +ysin θj

• Round ρ to the closest cell value, ρq

• Increment A(i,q) if the θj results in ρq

• After this procedure, A(i,j)=P means that P points in
the xy-space lie on the line ρj=x cos θj +ysin θj

• Find line candiates where A(i,j) is above a suitable
threshold

I have written the code to increment the indices of hough accumulator array.

My Code:

%%
% Hough Transform

[rows,cols] = size(E); % size of edged image

rhomax = ceil (sqrt(rows^2 + cols^2)); % maximum value of rho

rhomax = rhomax * 2; % multuplied with 2 to cater both +rho and -rho values

accu = zeros(rhomax, 101); % initialization of accumulaor array with rhomax and 101 possible theta values

for i =1:1:rows % for all rows

    for j =1:1:cols % for all columns

        if(E(i,j) ~= 0) % if the pixel is ONE (part of the edge) then process it

            for theta = 85.0:0.1:95.0 % for the theta values 85.0 - 95.0

                % computing the theta index for hough accumulator
                theta_index = 0;
                for k =85.0:0.1:theta
                    theta_index = theta_index + 1;
                end
                %---

                rho = round (i*cos(theta) + j*sin(theta)); % computing the value of rho

                % computing the rho index for hough accumulator
                rho_index = 1;
                if(rho < 0)
                    rho_index = rho_index + abs(rho);
                elseif(rho > 0)
                    rho_index = rho_index + (rhomax/2) + rho;
                elseif(rho == 0)
                    rho_index = rho_index + (rhomax/2);
                end
                %---

                accu(rho_index,theta_index) = accu(rho_index,theta_index)+1; % incrementing the corresponding index in accumulator

            end

        end

    end

end

Now I want to proceed further with the threshold value. I will be very grateful if someone tell me how should I select my threshold value ?

Well, it really depends on what you expect to find in the image.  The matlab function uses an example threshold of half the largest peak.

What sort of images are you using this on?

So, based on just finding vertical lines in your image, let's make some assumptions.

Your rectangle's minimum vertical side length is $L_{\rm min}$ pixels.
Your edge detection only finds 80% of each edge.

Then the minimum threshold will be
$$
0.8 L_{\rm min}
$$

You many need to modify this if, for example, your edges tend to be thicker or thinner than the 80% figure.


Classic Signal Processing vs Deep Learning / Machine Learning (DNN / ML) Based Signal Processing Are classic signal processing/statistics based approaches to optimum detection/estimation still relevant/important compared to ML based approaches using DL?
There was a time when speech processing was done using HMM. While not completely unused most speech processing now rests firmly in the DL area. I'm wondering if every area of classic signal processing will be taken over by these approaches over the classic feature finding historically used.

This is an open discussion.
In my opinion, any field we couldn't figure near optimal solution (And in most cases we only have Linear Optimal solution) using Deep Learning will replace classic methods given enough data to work with.
The power of data driven features vs. intuitive (Though mostly right yet never cover everything) will prevail.
I don't see it as a problem, Deep Learning isn't throw it all and it work machine, it still requires domain knowledge to effectively build the training process, the architecture and mostly the model to optimize. Which, if out think about it, just like in a classic method: Build test case, build a model, build a solver.
So one can see it as an evolutionary step where the optimization and the model are greatly improved yet the method has the same rational.
I believe near the hardware Signal Processing will stay as it is for longer time (Decimation, Filtering, etc...). Though in image processing DL is doing steps towards this as well. But in RF I think it will stay longer. But the next steps (For instance, Matched Filter) can be replaced (In Low SNR?) by DL.


Bio-signal discontinuity after FIR filtering I'm filtering respiration signal. I have a signal of duration equal to 10 minutes, so I decided to divided it into period of 10 s. I filtered each epoch with FIR filter, but I have discontinuities each 10s.
I attached an example.
I can not filter all the signal at once by the theory, so how can avoid discontinuity?

You need to maintain the filter state. For an FIR filter that's fairly easy since the state are just the last $N-1$ samples if your previous input, where $N$ is the FIR filter length.
Most filter functions have a built in option for doing this. If you don't have it just do the following:

Prepend you current frame with the $N-1$ last samples of the previous frame
Run the filter
Throw away the first $N-1$ samples of the output.


Equalization and then Sharpening vs Sharpening and then Equalization This was a question from an exam of digital image processing:

"Working with low contrast images, an elective method to enhance details and improve contrast is to apply in cascade an equalization of the histogram and a sharpening filter technique.
a) Does the order of application of the methods matter? Justify
b) Which method would you use first? Justify"

Could anyone help me? I can't find a clear answer in the textbooks.

Since the Histogram Equalization isn't Linear Operator the order does matter and the result of the doing one before the other changes the final result.
Sharpening would work best when all the Dynamic Range is used.
Hence I'd first maximize the Dynamic Range by Histogram Equalization and then would apply the sharpening.
You might also think of the noise when making this choice.


Calculating $E_b/N_0$ from given SNR ratio I have the following formula:
$$\frac{S}{N}=\frac{E_b \cdot R_b}{N_0 \cdot B }$$
Where $S/N$ is of course SNR in dB, $E_b$ is energy per information bit, $ R_b$ is information bit rate in [bit/s], $N_0$ is noise spectral density and $B$ is bandwitch in [Hz].
So let's say the $S/N$ ratio = $70\ \rm dB$

$R_b = 250\ \rm kb/s$
$B = 1\ \rm MHz$

So:
$$70\ \mathrm{dB}=\frac{E_b \cdot 250 000\ \mathrm{b/s}}{N_0 \cdot 1000000\ \mathrm{MHz}}$$
$$70\ \mathrm{dB}=\frac{E_b \cdot 1\ \mathrm{b/s}}{N_0 \cdot 4\ \mathrm{MHz}}$$
Now, of course, I can't simply multiply because this 4 is not on a logarithmic scale. But can I just logarithmize 4 like that? Won't there be a problem with units? How to do it to get $E_b/N_0$ ratio in decibels, because I keep doing something wrong.

HINT:
$$
\frac{S}{N}=\frac{E_b \cdot R_b}{N_0 \cdot B }
$$
With all units linear:
$$\implies \frac{E_b}{N_0} = \frac SN\cdot \frac BR_b\tag{linear}
$$
Now with $S/N$ in $[\rm dB]$ and $R_b$ and $B$ linear, you get $E_b/N_0$ in $[\rm dB]$ as follows:
$$\implies \frac{E_b}{N_0} = \frac SN +10\log_{10}\left(\frac BR_b\right)\tag{in [dB]}$$


How to compare the qualities of two interpolation (image resizing) algorithms? I am implementing image resizing algorithms (Bilinear, Bicubic, Lanczos and a few others).

How do I quantitatively compare them?

I am thinking of considering a large sample of images and running these algorithms on all of them and compare pair wise wrt a standard interpolation implementation. What metric shall I use here? Any other ways of doing the same?

Usually to test interpolation algorithm (for resizing, but also rotation and all sorts of "reversible" image transforms, you would apply the forward algorithm several times (ie resize to twice the size, take the result and resize it to 1.5x the size, etc .) and then apply all the backward transformations.

The more transforms you apply the more errors you will see and the more sensitive your comparison is. It can also represent your use case better (people pick an already resized image, resize it, and publish it)

Then you compute the max difference (or the standard deviation, or just visually compare) between the result and the original image.

In your case you can also apply all the forward transforms and then resize it once to the original size.


When and how does one use a Wiener filter? I'm trying to get some insight in this topic. As far as I understand, a determined signal enters a Wiener filter and the output is an estimate of some desired signal. Then, one can substract the desired signal to the output of the filter and calculate the estimation error. This diagram would represent what I just described above, where $x(n) = \hat{s}(n)$, the estimate of the desired signal $s(n)$, and $w(n)$ is some signal that has some correlation with $s(n)$:

I don't understand why I would try to estimate $s(n)$ if I already have it (I wouldn't be able to calculate the error $e(n)$ if I didn't have the desired signal).

The next diagram makes a bit more sense to me:

It would be a standard noise-reduction filter. A noisy signal comes in, a less-noisy one comes out.

There is a third case I found:

Here, one estimates the noise $v(n)$ to subtract it from a noisy signal $s(n)+v(n)$ and get a cleaner version of it, $\hat{s}(n)$. In this case, I have the same question as in the first one: why would I estimate the noise to subtract it from $s(n)+v(n)$ if I already have to know what the noise signal is in order to put it at the input of the filter?

So, in summary, I want to know if all of these cases are of use, and if they are equivalent in some sense. Also, I want to understand why they always estimate a signal that is already known, or if they don't do that and I'm not thinking correctly.

Perhaps a motivating radar/audio example would be an adaptive sidelobe canceller and adaptive noise cancellers. Rather than just show equations, let's walk through some text descriptions:

Let's say you are trying to record someone singing, but you're in a large auditorium. Your microphone is picking up the singer, but it's also picking up a lot of the reverberation. Now if we can characterize the reverberation somehow, we'll be able to remove it from our signal of the singer+reverberation. Enter another microphone (or a group of microphones), set up to just record the reverberation. With this setup, we can construct an adaptive filter that will remove the unwanted reverberation.

Another example would be some basic electronic protection for a radar system, where we want to remove jamming interference. In a radar system, you have a main channel which extracts the signal you're interested in. Sometimes, the case of some electromagnetic interference that originates from outside of the main antenna beam, i.e. from an antenna sidelobe, arises. This severely corrupts our main channel with noise, and may obscure target returns.

Intuitively, if we can somehow characterize just the interference, we can construct an adaptive filter that will remove the noise from the main channel. We can do this using single elements with isotopic (or nearly isotopic) responses. These auxiliary channels as they're called have very low gain compared to the entire antenna, and thus will receive only the interference and not a target which may be in the main beam.

The downside with wiener filters is that such a solution requires the noise to be stationary. Fortunately, we can use filter banks, time segmentation, and other signal processing techniques to make that possible.

Perhaps someone with a communication systems background could chime in and talk about adaptive equalizers, that's another application where you're trying to characterize some noisy channel, and you transmit a training signal, which the system knows about. The idea then is basically your first example: you know what you sent out, and you get some noisy version back. Using an adaptive filter, you can characterize that noise and now be able to remove it for anything else you send out (assuming the noise is perfectly stationary).

Hopefully that gives you some insight and motivation!


Fourier Transform Time Scaling I edited the initial question to make it general and thus applicable in more situations.

If the Fourier Transform of $x(t^2)$ is some function of $\omega^2$. Then what is the Fourier Transform of $x(bt^2)$, for $b>0$?

Is it $\dfrac{1}{b}$ times a function of $\left(\dfrac{\omega}{b}\right)^2$, or $\left(\dfrac{1}{\sqrt{b}}\right)$ times a function of $\left( \dfrac{\omega}{\sqrt{b}} \right)^2$, or neither?

Thank you for your time.

Let $y(t) = x(t^2)$, and $a$ be such that:

$$
\begin{align}
  y(at)
   &= x((at)^2) \\
   &= x(a^2t^2) \\
   &= x(bt^2)
\end{align}
$$
In other words, $a^2 = b$.

The Fourier Transform time-scaling property states that given the Fourier Transform pair:
$$
\begin{align}
  y(t) \leftrightarrow Y(w)
\end{align}
$$
and $a>0$, we have the relation:
$$
\begin{align}
  y(at) &\leftrightarrow \frac{1}{a}Y(\frac{w}{a})
\end{align}
$$

Thus,
$$
\begin{align}
  x(bt^2) &\leftrightarrow \frac{1}{a}Y(\frac{w}{a}) \\
          &\leftrightarrow \frac{1}{\sqrt{b}}Y(\frac{w}{\sqrt{b}})
\end{align}
$$

Now, if we define
$$
\begin{align}
Y(w) = G(w^2)
\end{align}
$$
Then,
$$
\begin{align}
Y(w/a) &= G(w^2/a^2) \\
Y(\frac{w}{\sqrt{b}}) &= G(\frac{w^2}{b})
\end{align}
$$

Yielding,
$$
\begin{align}
x(bt^2) &\leftrightarrow \frac{1}{a}Y(\frac{w}{a}) \\
        &\leftrightarrow \frac{1}{\sqrt{b}}G(\frac{w^2}{b}) \\
        &\leftrightarrow \frac{1}{\sqrt{b}}G((\frac{w}{\sqrt{b}})^2)
\end{align}
$$
Or, in prose (1/sqrt(b)) times a function of (w/sqrt(b))^2.


Envelope detector for non-zero mean signal (Matlab) I used hilbert transform to get an envelope signal.

envelope = abs(hilbert(inputSignal))

This works for zero-mean signal:

But, for nonzero-mean case, the output is not what expected.

If I do not want to remove the mean level of signal,
How can I extract the same envelope as in the first case?

You can calculate the mean of your nonzero signal, then substract it to the non-zero mean signal and calculate the envelope as you did before.

That will give you the envelope of the signal for the zero mean signal, but if you want to recover the envelope for the non-zero mean signal, then you just have to add to this envelope the initially calculated mean.


How to make sinewives with array input, for-loop functions and operate array with others [matlab]? Thanks to previous question, I got some clues to resolve this one.
But, I still get some problems here.

Ts = .001;
N = 1000;
t = (0:999)*Ts;
 make_an_array = [2 10 15 20 30 40 50 60 80 90 100 150 200 300 400];
 f = make_an_array;
for i=1:length(t)
input(i)=sin(2*pi*f(1)*t(i));
end
[b,a] = butter(4,.02,'High');
[b1,a1] = butter(2,0.12);
input = filter(b1,a1,input);
output = filter(b,a,input);
for i=1:length(output)
G(i) = 20*log(output/f(1));
end
semilogx(f,G,'x')
xlabel('frequency in Hz')
ylabel('Gain dB')
A = max(output);
grid on
disp(A)

This is a brief explaination about what I am trying to find:

1. Amplitude of the output with max operator (variable A)
2. Sine wives using input frequencies in an array (variable f) and use a for-loop (input (i))
3. Plot the 20 log of the output values against the frequency array using the semilogx.
4. I dont know about how to 'store the maximum values of the flter’s output in an array for plotting'.
5. And then, I got error message "Subscripted assignment dimension mismatch."
6. Part of codes is already on the book:
   [b,a] = butter(4,.02,'High');
   [b1,a1] = butter(2,0.12);
   input = filter(b1,a1,input);
   output = filter(b,a,input);

I think I got wrong on trying to find the value of G, how to operate value with array function (like in input(i) and G(i)).

I really don't like doing someone's homework so I'll drop the loops and make a few other modifications and leave it up to you figure it out.

I don't really understand why you filtered the input but I kept it in:

clc
clear all
close all
%%
Ts = .001;
N = 1000;
t = (0:999)*Ts;
%f = [2 10 15 20 30 40 50 60 80 90 100 150 200 300 400];
% f = make_an_array;
f=logspace(log10(2),log10(400),50);
%%
input=exp(1j*2*pi*f'*t);
[b,a] = butter(4,.02,'High');
[b1,a1] = butter(2,0.12);
input = filter(b1,a1,input);
output = filter(b,a,input);
G=output./input;
figure(1)
semilogx(f,20*log10(max(abs(G'))))
xlabel('frequency in Hz')
ylabel('Gain dB')
%A = max(output);
grid on
%disp(A)


Optimal Chebyshev FIR filters I have many information i want to be clear very well .

The defining characteristic of FIR filters optimal in the Chebyshev sense is that they minimize the maximum frequency-response error-magnitude over the frequency axis. In other terms, an optimal Chebyshev FIR filter is optimal in the minimax sense: The filter coefficients are chosen to minimize the worst-case error (maximum weighted error-magnitude ripple) over all frequencies.In another word, minimizing the worst case squared error induces a minimum Chebyshev error problem in some formulations.

In order to achieve an optimal equiripple complex FIR filter in the complex Chebyshev sense ,we can use an evolutionary algorithm.

What is the relation which allowed using the complex Chebyshev sense for optimisation equiripple complex FIR filter? What is effective evolutionary algorithm in this case?

The Chebyshev criterion minimizes the maximum weighted absolute error:

$$\epsilon=\max_{\omega\in\Omega}W(\omega)|H(\omega)-D(\omega)|\tag{1}$$

where $\Omega$ is the union of all frequency bands of interest, $W(\omega)$ is a positive weight function, $H(\omega)$ is the filter's complex frequency response, and $D(\omega)$ is the complex desired frequency response.

Note that in general, $(1)$ is a complex approximation problem, because $H(\omega)$ and $D(\omega)$ are generally complex-valued. However, in the special case of linear phase FIR filters, $(1)$ can be rewritten as a real approximation problem, for which a very efficient algorithm exists (the Remez exchange algorithm used in the Parks-McClellan method). In the general case of non-linear phase FIR filters or IIR filters, $(1)$ can be solved by (much) more complex general optimization algorithms. For FIR filters, the problem is still a lot easier than for IIR filters, and several algorithms have been suggested (cf. the thesis linked here and the references therein). One of the algorithms is implemented in Matlab's signal processing toolbox: cfirpm.m


I need help interpreting FFT diagram I need help with analyzing the computed FFT of my signal. I don't wanna bother by going into the details of how my signal is formed so I will just state which frequencies I know my signal has. I know that the fundamental frequencies of my signals are 1,325 Hz and 1,079 Hz.
After computing the FFT of my signal I got the following FFT diagram.

My first observation is that the fundamental frequencies I expect are not identical to the frequencies I calculated. Is that to be expected when doing FFT analysis of any signal? My expected frequencies are mentioned above, the similar frequencies to those are 1.28 ( close to 1,325 Hz) and 1.038 (close to 1,079)

Secondly, I thought that the fundamental frequencies in my signal should have the highest amplitudes from all other peaks. That is the case with the peak at 1.28 Hz but not the case with the peak at 1.038. What could be the reason for that?

Finally, I know that the other peaks represent the rest of the harmonics of my signal but I am not sure how to interpret them exactly? what do they exactly say about my signal?

Is that to be expected when doing FFT analysis of any signal?

No. Either your calculation or the implementation of your analysis is off.

Secondly, I thought that the fundamental frequencies in my signal should have the highest amplitudes from all other peaks.

Why did you think that? There are many signals where the fundamental is lower than the harmonics or even entirely absent.

Finally, I know that the other peaks represent the rest of the harmonics of my signal

These are not all harmonics of your expected frequencies. You have also peaks at 0.25Hz, 0.5Hz. Could also be noise spurs, intermodulation, etc.

what do they exactly say about my signal?

Without knowing what that signal actual is, it's impossible to even  guess.


Given a signal and its Fourier transform, find FS coefficient of the shifted sum of the signal Given $x_1(t),X_1(j\omega), x_2(t)=\sum_{k=-\infty}^{\infty}x_1(t-6k)$, find Fourier series coefficient of $x_2(t)$.

Looking up the FT table, I got $X_2(j\omega)=\sum_{k=-\infty}^{\infty}e^{-j\omega 6k}X_1(j\omega)$. Fourier transform can be represented as a summation of FS terms, so FS coefficient of $x_2(t)$, $a_k$, is $X_1(j\omega)$, is this correct?

Figure out that $x_2(t)$ is basically, sum of shifted copies of $x_1(t)$ which can be written as follows :
$$x_2(t) = x_1(t) * \sum^{\infty}_{k=-\infty} \delta(t-6k)$$, where, $*$ represents convolution operation. therefore, the Fourier Representation of $x_2(t)$ will be product of the fourier representation of $x_1(t)$ and the pulse train $\sum^{\infty}_{k=-\infty} \delta(t-6k)$.

Fourier Transform of a pulse train $\sum^{\infty}_{k=-\infty} \delta(t-6k)$ is given by the following :

$$\mathscr F \Big\{ \sum^{\infty}_{k=-\infty} \delta(t-6k) \Big\} = \frac{1}{6} \sum^{\infty}_{k=-\infty} \delta(f-\tfrac{k}{6})$$

Which basically means that Fourier Representation of $x_2(t)$ becomes :
$$X_2(f) = X_1(f) \cdot \frac{1}{6} \sum^{\infty}_{k=-\infty} \delta(f-\tfrac{k}{6})$$

So, $f = \frac{1}{6}$ is the fundamental frequency, since you have periodized $x_1(t)$ by $T = 6$. And, therefore, at all multiples of $\frac{1}{6}$, you will get a Fourier coeff, which will be equal to $\frac{1}{6} X_1(\frac{k}{6})$.

What this means is, $X_2(f)$ is nothing but sampled version of $X_1(f)$ at $f = \frac{k}{6}$ and hence, $X_2(f)$ is discrete. Which is evident from the fact that $x_2(t)$ was periodized in the first place.


How is LMS/FXLMS noise cancelling different than simple polarity inversion? Consider a noise-cancelling headphone:

If I have a noise signal from the outside world, mic it, flip the polarity by multiplying by -1 (with an op-amp or digitally), delay by the appropriate amount, and playback out of the speaker -- the summation of the noise and anti-noise at the ears should approach zero (with better performance at low frequencies).

How does FXLMS or LMS improve upon this? Why is polarity inversion ('phase cancellation') not as effective at cancelling noise, at least in the literature (e.g. Kuo 1994)?

Your polarity inversion method simply applies the same phase shift (delay) and gain to all frequencies, and is presumably not adaptive (or even closed loop). It may work reasonably well for some limited cases, where the phase and gain have been tailored for that particular environment, but any change to the environment will tend to reduce its effectiveness.

LMS however is an adaptive method which constantly "tunes" a filter to give you optimum cancellation using an appropriate gain and phase response which varies with frequency, and which adapts to a changing environment.


High Transmission Power : BPSK or QPSK? In BPSK,
the constellation consists of 2 points (equiprobable) spaced at distance of $\sqrt{E_b}$ each from the origin.

So the average power from the constellation can be obtained by:

For BPSK, $P = \frac{E_b+E_b}{2} = E_b.$
For QPSK, 4 points which are equiprobable are spaced at distance of $\sqrt{2E_b}$ from origin each, that is they are placed on vertices of a square of side $2\sqrt{E_b}$. So, for QPSK, $P = \frac{2\cdot4\cdot E_b}{4}=2E_b.$

So, does this mean that the power required for QPSK is more than BPSK when the probability of error is same for both?

Note: $E_b$ is the energy per bit and $P$ is the power which is assumed as proportional to the distance of the signal point from origin in the constellation.

That's right. When we compare BPSK and QPSK at equal energy per bit, we find that QPSK requires twice the power of BPSK. From wikipedia:

  Although QPSK can be viewed as a quaternary modulation, it is easier
  to see it as two independently modulated quadrature carriers. With
  this interpretation, the even (or odd) bits are used to modulate the
  in-phase component of the carrier, while the odd (or even) bits are
  used to modulate the quadrature-phase component of the carrier. BPSK
  is used on both carriers and they can be independently demodulated.

  As a result, the probability of bit-error for QPSK is the same as for
  BPSK:

  $P_b = Q\left(\sqrt{\frac{2E_b}{N_0}}\right).$

  However, in order to achieve the same bit-error probability as BPSK, QPSK uses twice the
  power (since two bits are transmitted simultaneously).


What is the antialiasing algorithm employed by video cards? What is the algorithm employed by video cards when one talks about for example 8xAA?
I thought the algorithm was about looking at neighboring pixels for correction. Yet, I hear that the implementation is just upsampling followed by downsampling.

There are various algorithms to do anti-aliasing. I think that the one that you mean is called super-sampling. The idea is that the scene is rendered in x8 resolution of what you eventually see. Then it is low-pass filtered, to eliminate aliasing of higher frequencies in the next step, which is down-sampling x8.

This is described in detail here.

There is a more general discussion of AA algorithms here, and another type of AA algorithm called Multisample AA here.


Fast fourier transform and frequency domain spike in wrong area I used MATLAB to produce a simple sine wave over time:

frequency=7
sampling_frequency=1000;
total_time=1-1/sampling_frequency;
step=1/1000;
theta=0;
time=0:step:total_time;

y=Amplitude*sin(2*pi*frequency*time+theta)

I then used the following code to fast fourier transform:

Y=fft(y);
ploty=Y/length(y);
stem(abs(ploty))

I was expecting to see one spike at 7 instead I am getting 2 and not even at 7.
One is at 8 and the other at 994. I am totally confused as to why.

You see one at index 8 because MATLAB's indexing is one-based, not zero-based as the equations that describe the DFT are. You get two spikes because your input signal is real and is therefore conjugate-symmetric in the frequency domain. If you work out the math, you'll see that the Fourier transform of a $\sin$ has two spikes symmetric about zero frequency. In order to get a single spike, you need to use a complex exponential function at the same frequency instead.


does ifft of higher order equal upsampling signal I have a signal with length 2000, and am taking an fft by 2048, after I perform analysis and some operation on signal in frequency domain I tried to convert it back to time domain. I received a signal length 2048, looking similar to what I would expect but not the same. I know that taking higher point fft simply generates more samples of the DTFT we would expect from the original signal, it does not add new information. However, when I take ifft my data length increases, I am assuming it must time domain interpolation between signal and output is same as if I would upsample my expected 1000 point signal (the operations I perform can be ignored, in my case if this is true for simply taking higher fft and then ifft of signal as is, then its true for me). I am not sure if this assumption is true, as I have yet to grasp how the IDFT works when length in frequency domain changes. I am performing analysis using MATLAB and am aware of the sampling rate and time. Can anyone please share some insight on if my assumption is true and any mathematical insight into it or point me to some freely available resource?

The longer FFT simply pads zeros to the original time domain signal. The inverse FFT should then be the original signal with zeros appended to the end.
The OP is correct that taking a longer FFT interpolates samples in frequency, and specifically provides more samples of the Discrete Time Fourier Transform. Why this occurs is clearer by observing the difference between the DFT and DTFT formulas as given below:
DFT:
$$F(k) = \sum_{n=0}^{N-1}x[n]e^{j k \omega_nn }$$
Where $\omega_n = 2\pi /N$
DTFT:
$$F(\omega) = \sum_{n=-\infty}^{\infty}x[n]e^{j \omega n }$$
For $\omega$ continuous from $0$ to $2\pi$ (for unique span)
Notice the two key differences: The time increment n in the DTFT extends to plus and minus infinity, and the frequency domain variable $\omega$ is a continuous function. As we add zeros to x[n] in the DFT, we approach the DTFT and ultimately for the same $x[n]$ all results of the DFT as more and more zeros are added will be samples on the continuous DTFT function.
In general related to this it is worthwhile to understand a universal property of all the Fourier Transforms related to signals that repeat or not (periodic or aperiodic) and continuous or sampled signals:
If a signal is continuous in one domain (time or frequency) it must be aperiodic in the other domain. And the counter-case is also true; if the non-zero values of a signal are discrete in one domain, it must be periodic in the other if we allow the domain to extend toward infinity. We see that immediately with the premise of the Fourier Series Expansion in that any single-valued analytic waveform over a finite interval can be represented by a discrete and infinite number of frequencies; adding together all those frequencies in the time domain over an infinite time duration will result in a periodic waveform.
With that observe the Continuous Time Fourier Transform (CTFT) as shown in the slide below. Regardless if the waveform itself is zero valued over a portion of the interval, the key point is the Transform itself uses a time axis that extends from minus infinity to plus infinity, and for that reason the frequency domain MUST be a continuous function: Frequency resolution is the inverse of observation time. If we have infinite time to observe a stationary signal, then we will have infinite frequency resolution. So the CTFT is of a signal in time that is continuous and aperiodic. Therefore the result in frequency is aperiodic and continuous.

Now compare this to the Fourier Series Expansion (FSE) as shown in the slide below, of a time slice of the same time domain signal. As previously explained, if the time is extended beyond the time interval from $0$ to $T$ the resulting waveform will be periodic. The FSE is specifically done over the interval from $0$ to $T$, but periodicity is implied as the exact same result in frequency would occur at the non-zero values if time was extended to minus and plus infinity. The frequency axis for the FSE is discrete since the results are defined for $n \omega_o$ with integer n. What also would occur if time was extended to infinity as explained for the CTFT is the frequency axis would become a continuous function, but the new interpolated values between the original samples would all be zero valued. The FSE of a signal in time from 0 to T is the same as a signal that is periodic in time over that same period. The signal in time is continuous and periodic, therefore the result in frequency is aperiodic and discrete (non-zero values are discrete).

Now compare this to the Discrete Time Fourier Transform (DTFT) as in the slide below for the same time domain waveform.  The only difference from CTFT is that the signal is now sampled in time. Thus the signal that is now discrete and aperiodic in time will be periodic and continuous in frequency: The same frequency result in the CTFT becomes periodic in the DTFT.

Finally the Discrete Fourier Transform (DFT) combines both effects of the FSE and DTFT; time is discrete as in the DTFT and restricted to a finite interval as in the FSE, with implied periodicity in that we can equally extend time to plus or minus infinity without changing the non-zero results in frequency. Thus the DFT for the same time domain waveform as the cases above would appear as in the slide below. The signal that is discrete and periodic in time would be periodic and discrete in frequency.


Can I create my IF signal on SDR before bringing data to PC? I'm relatively new to SDR's and I'm trying to determine what sort of BW I need for my laptop to get the best resolution possible for an FMCW system. I have an X310 Ettus with a UBX160 daughtercard which should give about 1 meter of resolution theoretically. Anywhere close to this will probably work just fine for my application. But at the moment I only have a PC capable of a 1 Gig ethernet hookup. The fact is though this will be a pretty short range system, and I expect a pretty narrow range of results (approx bandwidth of a few hundred kHz). Is there a way to mix the echo with the FMCW waveform on board the ettus and just deliver the sub MHz band to the laptop? Or do I really need to shell out all that cash to get the range resolution this unit could deliver?
Sorry this is a newb question, I've been playing around with GNU Radio and uhd a little and I'm leaning towards expecting a "no" on this question - but would love some input from the pros out there. Thanks in advance!

The theoretical data rate to transmit a signal with bandwidth $W$ Hz and resolution $R$ b/sample is simply $2WR$ b/s. In practice, you will need somewhere between that and $4WR$.
One way to reduce the required rate is to preprocess the signal in the SDR's FPGA. This will require delving into Verilog, integrating your hardware with what is already there, and very likely simulating and debugging it. What can be done also depends on how much space is left on the particular FPGA in your device.


How can I make a mmWave frequency selective channel in matlab? $$H_d=\sum_{\mathcal{\ell}}^{N_p}\alpha_\ell p(dT_s-\tau_\ell)\mathbb{a}_R(\phi_\ell)\mathbb{a}^*_T(\theta_\ell) $$

I'm trying for simulation of a frequency selective channel in matlab. my problem is the raised cosine filter. without the raised cosine filter the matlab code is as follows. but in presence of pulse shaping filter I confused. can any body help me?

EDITED:

I've done all of Mr. Carlos Danger advises and this is my final code. but I got some errors. can any body help me?
 this is my RRC code:

function [res]=Root_Raised_Cosine(Sample_Rate,Roll_Off_Factor,time)

Sample_Interval = 1/Sample_Rate;
a = Roll_Off_Factor;
index = [-30*Sample_Interval:Sample_Interval:30*Sample_Interval];
%This can be increased or decreased according to the requirement

p = zeros(1,length(index));
for i = 1:1:length(index)
    if index(i) == 0
        p(i) = Sample_Rate*((1-a)+4*a/pi);
    else if index(i) == Sample_Interval/(4*a) || index(i)==-Sample_Interval/(4*a)
           p(i) = Sample_Rate*a/sqrt(2)*((1+2/pi)*sin(pi/(4*a))+(1-2/pi)*cos(pi/(4*a)));
          else
            p(i) = Sample_Rate*((sin(pi*index(i)*(1-a))+4*a*index(i)...
                .*cos(pi*index(i)*(1+a)))./(pi*index(i).*(1-(4*a*index(i)).^2)));
         end
    end
end
response = p./sqrt(sum(p.^2)); %Normalization to unit energy
t = find(abs(index-time)== Sample_Interval) ;
res = response(t);
end

and this is my channel.

clc; clear all;

%--------------------------System Parameters------------------------------%

Num_BS_Antennas = 16; % BS antennas
Num_BS_RFchains = 4;  % BS RF chains
BS_Antennas_Index = 0:1:Num_BS_Antennas-1; % Indices of the BS Antennas

Num_MS_Antennas = 16; % MS antennas
Num_MS_RFchains = 4;  %MS RF chains
MS_Antennas_Index = 0:1:Num_MS_Antennas-1; % Indices of the MS Antennas

Num_Qbits = 7; % No. of Quantization bits
Step = 2^Num_Qbits;
Step_Index = 0:1:Step-1;
PS_Quantized_Angle = Step_Index*(2*pi)/Step; % Quantized angles in PSs

%---------------------- Hybrid Precoders and Combiner---------------------%

for i = 1 : 1 : Num_BS_Antennas
for j = 1 : 1 : Num_BS_RFchains
    uniform_index = ceil(length(PS_Quantized_Angle)*rand);
    F(i,j) =
sqrt(1/Num_BS_Antennas)*exp(1j*PS_Quantized_Angle(uniform_index)); % Hybrid
Precoder
end
end
clear i j
for i = 1 : 1 : Num_MS_Antennas
for j = 1 : 1 : Num_MS_RFchains
    uniform_index = ceil(length(PS_Quantized_Angle)*rand);
    W(i,j) =
sqrt(1/Num_MS_Antennas)*exp(1j*PS_Quantized_Angle(uniform_index)); % Hybrid
Combiner
end
end

%--------------------------Channel Parameters ----------------------------%

Delay_Tap_Length = 4;
Sample_Rate = 1760*10^6; % Sampling rate set to be 1760 MHz
Sample_Interval = 1/Sample_Rate;
delay = 0 : Sample_Interval : (Delay_Tap_Length-1)*Sample_Interval;
Num_paths=3; % Number of channel paths
Roll_Off_Factor = .8;
% Channel parameters (angles of arrival and departure and path gains)
AoD=rand(1,Num_paths); %Sin of AoD in the range og [0,1)
AoA=rand(1,Num_paths); %Sin of AoA in the range og [0,1)
alpha=(sqrt(1/2)*sqrt(1/Num_paths)*
(randn(1,Num_paths)+1j*randn(1,Num_paths)));

% Channel construction
 d = 3;
Channel=zeros(Num_MS_Antennas,Num_BS_Antennas);
for l=1:1:Num_paths
Abh(:,l)=sqrt(1/Num_BS_Antennas)*exp(1j*pi*BS_Antennas_Index*AoD(l));
Amh(:,l)=sqrt(1/Num_MS_Antennas)*exp(1j*pi*MS_Antennas_Index*AoA(l));
res = Root_Raised_Cosine(Sample_Rate,Roll_Off_Factor,d*Sample_Interval-
delay(l));
Channel=Channel+
  (sqrt(Num_BS_Antennas*Num_MS_Antennas)*alpha(l)*Amh(:,l)*Abh(:,l)');

end

Your code appears to be summing up amplitudes over all the paths to simulate Rayleigh (flat) fading.  The assumption of flat fading is that the delay spread of the channel is much less than the symbol period, so all the symbols add up roughly at the same time (1).

To make the leap to frequency-selective fading, you need to account for the pulse shape in your summation, i.e. you need to have some vector of sample indices n, a function p(n) that returns the pulse amplitudes at times n*Ts, and multipath delay ndelay. Something like this (pseudo-MATLAB code):

n = 0:maxn;
Npath = 4;
ndelay = [2, 5, 3, 7]; % length(ndelay) == Npath
sumval = 0;
for i = 1:Npath
    sumval = sumval + otherstuff*p(n - ndelay(i));
end

You must specify ndelay ahead of time manually or generate it according to some probability distribution.

References:
(1) Goldsmith, Wireless Communications, chapters 2-3, 2005.


FM synthesis fuzzy sound I'm trying to reproduce this !!!WARNING LOUDNESS!!! sound in C. In the website, it says to use Carrier: 500 Hz, modulator frequency: 1 Hz, modulator index: 500, however I'm getting lots of clicks and fuzzing (along with the right sound) with this code:

#include <stdio.h>
#include <math.h>
#include <string.h>

#define DURATION    20
#define RATE        44100
#define NUMCHANS 1

#define NUMSAMPLES  (DURATION*RATE)
#define NUMSAMPLESC (NUMSAMPLES*NUMCHANS)

#define BYTES_IN_HEADER (11*sizeof(int))
#define BYTES_IN_BUFFER (NUMSAMPLESC*sizeof(short) + BYTES_IN_HEADER)

#define DELTAT (1.0/RATE)

#define M_PI 3.14159265359

static const int wavHeader[] = {
0x46464952,
NUMSAMPLESC*2+36,
0x45564157,
0x20746D66,
16,
1|(NUMCHANS<<16),
RATE,
RATE*NUMCHANS*sizeof(short),
(NUMCHANS*sizeof(short))|((8*sizeof(short))<<16),
0x61746164,
NUMSAMPLES*sizeof(short)
};

static short buffer[BYTES_IN_BUFFER/sizeof(short)];

int main(void) {

 unsigned int i;
 double t = 0.0;
 FILE *fp;

 memcpy(buffer, wavHeader, BYTES_IN_HEADER);

for(i = BYTES_IN_HEADER/sizeof(short); i < BYTES_IN_BUFFER/sizeof(short); i++, t += DELTAT) {
    buffer[i] = 100 * sin(2*M_PI*500*t + 500*sin(t*2*M_PI*1));
}

 fp = fopen("test.wav", "w");
 fwrite(buffer, 1, BYTES_IN_BUFFER, fp);

fclose(fp);

return 0;
}

I have tried using single precision floats, reversing endianness, but with the same result. The only thing that seems to remove the clicks and fuzzing to a certain extent (but they're still there) is if i copy the lower 8 bits of each sample into the upper 8 bits.

My sound card is integrated from a P43DE3 motherboard and can handle audio CDs well, and my OS is Windows 8.

How do I reproduce that sound without noise?

EDIT

It seems to produce the correct output when I run it on a Mac. Perhaps and endianness error? Will investigate further.

OK, I figured out what the problem was using a hex editor. The compiler on Windows was seemingly introducing random 0x0D throughout the file. That's because of the mode I was opening the file with, "w". fwrite converted every instance of '\n' into "\r\n". All I had to do was change the file opening mode to "wb+", and it sounds exactly like the linked audio file, after changing the carrier amplitude to 2e5.


Does accurate translation estimation using a camera matrix mean the orientation estimation will also be accurate? Pose Estimation of single camera and planar pattern:

If I have an appropriate camera matrix with intrinsic parameters that allows me to calculate translation of the camera relative to a chessboard, for example, and it does so very accurately, does that mean my rotation calculations will also be accurate?

This may also be considered in the context of the camera calibration documentation for opencv.

References to your answers would also be appreciated.

EDIT:
thanks to @tbirdal for the links
"Pose and position" is an incorrect phrase. Pose is a combination of position and orientation.

If you pose estimation algorithm carefully minimizes the reprojection error using the standard non-linear techniques e.g., Levenberg Mardquardt, and your corner detection scheme is sufficiently accurate in the subpixel level, then yes. I would expect the residuals to be very small. Of course take into account the distortion parameters; they have to be calibrated as well. If you are using OpenCv and referring to the standard functions such as findExtrinsicCameraParam, then you are good to go.

The problem is also known as PnP due to the correspondences. Lepetit et. al. review it nicely here.

Also, Grest et.al provide an overview of 3 methods here.


Transform a Continuous Time Differential Equation of a Low Pass Filter to Discrete Time Assuming I have a linear differential equation of first order low pass filter:
$$RC*y(t)'+y(t) = x(t) $$
And Assuming i know now that the input x(t) is being sampled at a rate of   $T_s$

how can now switch the equation to describe in the discrete time form?
I know that the answer is:
$$x[n] = y[n] + RC/T_s* (y[n]-y[n-1])$$

but i can't figure how to get there.
I know that there is a strong relation between DTFT and fourier transform but i couldn't figure how to use it.
Thanks

The answer comes from a discrete approximation of the derivative operator.
Sometime it is called "Finite Differences".

Look at $ \frac{ y \left[ n {T}_{s} \right] - y \left[ n {T}_{s} - {T}_{s} \right] }{ {T}_{s} } $, what happens when $ {T}_{s} \to 0 $?
We got the definition of the drivative.

Now, since we sample at interval of $ {T}_{s} $ we can agree that $ y \left[ n {T}_{s} \right] = y \left[ n \right] $ and $ y \left[ n {T}_{s} - {T}_{s} \right] = y \left[ n - 1 \right] $.

Do that and you get:

$$ R C {y}' \left ( t \right ) + y \left ( y \right ) = x \left ( t \right ) \approx R C \frac{ y \left [ n \right ] - y \left [ n - 1 \right ] }{ {T}_{s} } + y \left [ n \right ] = x \left [ n \right ] $$

Now simple adjustment and you have what you need.


While converting discrete filter specs into continuous filter specs, why do we multiply gain by T As I know so far that when we trying to design a DF(discrete filter) with specific specs "in terms of $\omega$", we take these specs and map it into continuous domain $H(j\Omega)$ using a certain method as "Impulse invariance" and design it by continuous design methods like "Butterworth".

As impulse invariance method says:
$$\omega = \Omega.T$$
and that part is fine to map $\omega_p(passband),\omega_s(stopband)$

But do we have to map also the gain of the filter?

i.e. $\delta_1,\delta_2$ and multiply them by $T$

it is, in my biased opinion, because nearly every textbook with the sampling theorem puts the $T$ gain factor in the wrong place.

we've had nearly the same question asked a couple decades ago at the USENET comp.dsp newsgroup.  it was someone asking how to design a LPF reconstruction filter with a passband gain of $T$.  they were wondering if it made a difference if the gain was 20 (for microseconds) or 0.00002 (for seconds).

so Dilip here describes the sampling theorem making the same mistake (in judgment for where the $T$ should go). to fix this, i think that this old Wikipedia version of the Nyquist-Shannon Sampling theorem scales it correctly.

the reason is this, the sampling function should be the dirac comb scaled with $T$ so that the ideal brickwall reconstruction filter has a gain of 1 (or 0 dB) in the passband.  the reason for this is that the Fourier series coefficients of the properly-scaled dirac comb, all of those coefficients are 1.  not $\tfrac1T$

$$ T \sum\limits_{n=-\infty}^{\infty} \delta(t - nT) = \sum\limits_{k=-\infty}^{\infty} 1 \ e^{j (2 \pi k/T) t} $$

then to really understand how things are practically reconstructed, the sampled signal

$$\begin{align}
x_\text{s}(t) &= x(t) \, \cdot \ T \sum\limits_{n=-\infty}^{\infty} \delta(t - nT) \\
     &= \ T \sum\limits_{n=-\infty}^{\infty} x(t) \, \cdot \delta(t - nT) \\
     &= \ T \sum\limits_{n=-\infty}^{\infty} x(nT) \, \cdot \delta(t - nT) \\
\end{align}$$

but because this is also true (from the Fourier series of the scaled dirac comb):

$$\begin{align}
x_\text{s}(t) &= x(t) \, \cdot \ T \sum\limits_{n=-\infty}^{\infty} \delta(t - nT) \\
     &= x(t) \, \sum\limits_{k=-\infty}^{\infty} e^{j (2 \pi k/T) t} \\
     &= \sum\limits_{k=-\infty}^{\infty} x(t) \,  e^{j (2 \pi k/T) t} \\
\end{align}$$

and that means that the spectrum of the sampled signal is

$$\begin{align}
X_\text{s}(f) &= \sum\limits_{k=-\infty}^{\infty} X(f) \Bigg|_{f\leftarrow f+k/T} \\
     &= \sum\limits_{k=-\infty}^{\infty} X\left( f-\tfrac{k}{T} \right) \\
\end{align}$$

to recover $X(f)$ (which recovers $x(t)$), we eliminate (using an ideal brickwall filter) all of the terms except the $k=0$ term.

$$\begin{align}
X(f) = \operatorname{rect}(fT) X_\text{s}(f) &= \operatorname{rect}(fT)  \sum\limits_{k=-\infty}^{\infty} X\left( f - \tfrac{k}{T} \right) \\
&= \sum\limits_{k=-\infty}^{\infty} \operatorname{rect}(fT)   X\left( f - \tfrac{k}{T} \right) \\
\end{align}$$

in the last summation all terms are zero except the $k=0$ term which remains untouched.

the $ \operatorname{rect}() $ function is simply

$$ \operatorname{rect}(u) \triangleq \begin{cases}
1 \quad & \mbox{if } |u|<\tfrac12  \\
0       & \mbox{otherwise}
\end{cases} $$

then to understand how a Zero-order hold is modelled, you send the ideally dirac sampled signal $x_\text{s}(t)$ through a filter that has this for an impulse response:

$$ h(t) =  \frac{1}{T} \operatorname{rect}\left( \frac{t - \tfrac{T}2}{T} \right)
 = \begin{cases}
\frac{1}{T} & \mbox{if } 0 \le t < T  \\
0           & \mbox{otherwise}
\end{cases} $$

which results in this piecewise constant output

$$\begin{align}
y(t) &= h(t) \circledast x_\text{s}(t) \\
 &= \ T \sum\limits_{n=-\infty}^{\infty} x(nT) \, \cdot \frac{1}{T} \operatorname{rect}\left(\frac{t - nT - \tfrac{T}2}{T} \right) \\
 &= \sum\limits_{n=-\infty}^{\infty} x(nT) \cdot \operatorname{rect}\left(\frac{t - nT - \tfrac{T}2}{T} \right) \\
\end{align}$$

that shows the correct scaling for the Digital-to-Analog converter (DAC) output.  and the frequency response corresponding to the ZOH is

$$ H(f) \, = \mathscr{F} \{ h(t) \} \,= \frac{1 - e^{-j 2 \pi fT}}{j 2 \pi fT} = e^{-j \pi fT} \operatorname{sinc}(fT)  $$

and that is properly scaled so that the ZOH gain at DC is 0 dB with no dimensional constant like $T$.

so the reason for your question is that nearly every textbook does this wrong by leaving the $T$ factor out of the sampling function and putting it into the passband of the reconstruction filter.


What is the flow of data from a wifi receiver to the application layer? I am not sure if I should ask this here but I am researching to build a receiver similar to a wifi one. I don't know here to look to find the flow that the data goes through.

For example, what I believe right now is that

the data stream gets captured by the antenna.
Now you have data to work with and it either gets put into packets in the receiver or in the computer.
Once in the computer, I assume it goes through a device driver for the wifi (usb for example).
Then from there I don't know where it goes.

If you can point me to a
   resource that would help with creating and programming a receiver for
   a different spectrum, that would be awesome. Thanks

the data stream gets captured by the antenna.

No, that's not right. The antenna doesn't "capture data"; it just converts  electromagnetic waves to electrical signal that your receiver can work with. It's still analog RF oscillations, and not digital "data".

The receiver then goes and converts that RF to baseband, i.e. moves it from 2.4 GHz to 0 Hz, for example, resulting in analog complex baseband signal.

That is then simply sampled (like a microphone-generated electrical signal with a sound card, but thousands of times faster) at regular intervals: it now becomes a digital signal, i.e. an "infinite" sequence of complex numbers in a chip.

These numbers are then analyzed by algorithms that – all within your Wifi chip –

look for the presence of a Wifi transmission
recognize the start time and exact frequency of that
correct tat
divide the result in many parallel, smaller chunks of signal
on each of them, correct influences of the transmission over the air
then, decide which symbol (i.e. a complex number) was originally sent
convert the resulting symbol stream to a stream of code bits
processing the stream of bits with an error-correcting code
passing the resulting info bits on to the next layer

That next layer typically is still on the same IC; it now understands what bits belong to a data packet, checks whether the error correction has been successful (checksum), checks whether the packet is relevant to the receiver at all, and then hands it of to the host (over whatever link you use – USB, PCIe, SPI...).

There, a wifi driver takes the "surviving" packets and analyzes them to see whether they contain data that is of any relevance to the network stack of your computer.

  Then from there I don't know where it goes.

Network stack.

  If you can point me to a resource that would help with creating and programming a receiver for a different spectrum

You need to learn about Software-Defined Radio, then (book example). That will make little sense if you haven't had basics of communications; the book I linked to tries to bring these basics to you, but it's certainly not the only book. If you're currently enrolled at university, you'd want to study their course on "basics of communications technology" or so.

Your Wifi card will not help you very much - it's not meant to be anything than a wifi receiver.

So, really, the hard part is not so much letting the data flow – your network stack of the operating system you're using will take care of that – but coming up with a receiver, which in essence is nothing but an estimator for the data sent by someone far away, whose analog signals were subject to a fiercely antagonistic channel.


Wavelet transform 3D plot for CoP I'm trying to perform wavelet transform and make a 3D plot like :

With the wavelet transform function :

$$
\textrm{CWT}_x^\psi (\tau, s)=\frac{1}{\sqrt{\lvert s\rvert}}\int x(t)\psi\left(\frac{t-\tau}{s}\right)dt
$$

Where $t$ is translation and $s$ is scale.

These are MATLAB and Python functions for wavelet transform:

MATLAB:[coefs,sgram,frequencies] = cwt(x,scales,wname, samplingperiod,'scale')
Python:
pywt.wavedec(data, wavelet, mode='sym', level=None)
(cA, cD) = dwt(data, wavelet, mode='sym')
scipy.signal.cwt(data, wavelet, widths)

I know to analyze the signal I have to move the wavelet (translation) to cover all of the signal. The functions of both MATLAB and Python need scales as parameter but there is nothing about translation. The $x$-axis is scales, the $y$-axis is translation.

I assumed $z$ is 2D (surface) because I need the coloring but I dont know what it is. Is it coefficients ?
what are approximation and detail coefficients ?
And what's translation? Is it one to length of my data array (number of data points) ?

I'm new in DSP and I'm confused if anyone can help me I'll appreciate it.

Update :

my data:

0.01009
0.010222
0.010345
0.010465
0.010611
0.010768
0.01089
0.011049
0.011206
0.011329
0.011465
0.011613
0.011763
0.011888
0.012015
0.012154
0.012282
0.012408
0.012524
0.012664
0.012791
0.012918
0.013043
0.013157
0.013284
0.0134
0.013516
0.013666
0.013793
0.013909
0.014024
0.014143
0.014271
0.014398
0.014515
0.014618
0.014722
0.01484
0.014957
0.015075
0.015192
0.015298
0.01539
0.015493
0.015598
0.015695
0.015776
0.015884
0.015978
0.016073
0.016157
0.016254
0.016363
0.016473
0.016572
0.016694
0.016803
0.016913
0.017021
0.017154
0.017242
0.017342
0.01745
0.017555
0.017648
0.017743
0.017851
0.017957
0.018065
0.018194
0.01831
0.018439
0.018582
0.018713
0.018843
0.018995
0.019137
0.0193
0.019464
0.019625
0.019781
0.019945
0.020124
0.020304
0.020447
0.020619
0.020762
0.020931
0.021088
0.021254
0.021398
0.021531
0.021648
0.021814
0.021965
0.022109
0.022251
0.022408
0.022563
0.022748

I used morlet wavelet , 1:150 scale and I got this result:

I get trough at scales 50, 150 , 250 , ... and peaks at 100, 200, 300 , ...
Why?

Basically, I see a plot of a 2D function of discretized scale and translation parameters. Instead of a smooth 2D surface, it looks like 1D plots of coefficients  at all scales $s_n$, put behind each other along each location $\tau_n$ on the translation axis. And each 1D plot is colored in a level-set fashion: the  "vertical" coloring below the 1D curve is related to the amplitude, potentially with color cycling: red for low amplitudes, then yellow, green, cyan, blue, magenta, red (for high amplitudes). This is apparently an instance of a waterfall plot:

  curves are staggered both across the screen and vertically, with
  'nearer' curves masking the ones behind

with amplitude coloring. So:

Coefficients (absolute value) gives you the height of the 1D curve (top) and the coloring below.
You don't have approximations and details with CWT. This is not DWT. "Only" low-scale to high-scale "detail" coefficients. There is not father wavelet or scaling function in that case.
Yes, without special settings, standard wavelet codes compute coefficients at each sample.

Alternatively, you can draw ribbon plots, with a Matlab code for the color (ribboncoloredZ.m):


Precision of Computer Vision algorithms Let's say the task is to determine element position on image.
First very important thing is correct detection of object then some algorithms of calculating position are used (for examble blob analysis). Everything depends on multiple things (detection correctness, used algorithms etc.)

Lets assume we have callibrated image and know error given by callibration.
What are the methods to calculate reliably precision of computer (and machine) vision algorithms?
Can it be done analiticaly or only by experiments and tests?

The question addres cases when we detect element position and also other computer vision problems.

I want to get references to problems which are related to computer/machine vision especially element position detection and present some correctness computations either analyticaly or experimental approach to show this precission.

Also suggestions how to improve this question are welcomed.

For example, Hartley & Zisserman suggest using preconditioning prior to homography estimation, because taking direct matrix inverse can lead to huge errors or instabilities.
This applies to any numerical method working with matrix inverse.

Feature detection algorithms often use sub-pixel approximation of interest point location.

Most books discussing numerical methods also deal with their stability analysis.

Sometimes you need to do some statistics to analyze precision and accuracy of your estimator (be it a least-squares estimator or maximum likelihood estimator). This is useful in algorithms like RANSAC, which deal with outliers. You would also like to know, how well the estimated transform fit your data and possibly discard results that are too inaccurate.

When working with finite differencing or doing some filtering, a slight Gaussian blurring is done to remove noise, which would otherwise cause huge errors in second derivatives.

Some problems in computer vision are ill-posed. A method of regularization (such a Tikchonov regularization) is necessary in order to solve them. Examples where this is necessary include computing anisotropic diffusion.


Help in deriving the update rule for equalizer Considering an FIR system excited by an input $s_n$. If the input is derived from a mapping, and a +1/-1 BPSK signal is derived. Considering the Constant Modulus Algorithm (CMA). The channel is modeled as an FIR. If the received signal is fed to the equalizer, the output is : $y_n = w^T(n) x_n$ where $w^T$ is the equalizer tap vector. The update equation for weights $w_{n+1} = w  - \mu \Delta J(w)$

Will there be a modification in the derivation of the CMA and the update rule?

The mapping from $z_n$ to $s_n$ is not bijective, so it can't be inverted. This means that even in the transmitter, $z_n$ cannot be recovered from $s_n$. This is why the constant modulus algorithm (CMA) need not (and cannot) be modified to estimate $z_n$. It can only be used to estimate $s_n$, and for this purpose no modification is necessary.

The thesis you referred to in a comment is mainly about blind signal separation. The update rules (3.16) and (3.17) apply to the iterative maximization of the log-likelihood function in order to compute the maximum likelihood solution of a source separation problem. The nonlinear vector-valued function $\mathbf{g}(\mathbf{u})$ depends on the probability density functions of the different source signals. The underlying model and the problem formulation are different from the model and the problem underlying the CMA, so there is no direct relation between the two, and the update rules mentioned above are not applicable to the CMA.


Why does a wave continuously decreasing in frequency start increasing its frequency past the half of its length? I'm trying to programmatically generate a wave (sine or square) with a continuously decreasing frequency.
To do so, I simply multiply the starting frequency for a decreasing value, that linearly goes from 1 to 0, at intermediate steps throughout a given time span.
I was puzzled to see (and hear) that, exactly at the half of the time span, the frequency started to increase again.
I checked the frequency value but everything is ok.
I can use only the first 'half' of the wave, but it would be great to understand what I'm experiencing.
Any help appreciated !
Thanks

Make sure that your frequency doesn't reach values below 0 or above the half of your sample rate.

Please post more information/code about how you generate your waveform! Chances are you are not doing it correctly.

For example, if you want to generate a sine wave with a time-varying frequency $f(t)$ (for example to implement frequency modulation), generating something like:

$$y(t) = \sin (2 \pi t f(t))$$

is wrong, because your instantaneous frequency is:

$$f_i(t) = \frac{1}{2 \pi}\frac{d \phi(t)}{dt}$$

And:

$$\frac{1}{2 \pi}\frac{d 2 \pi t f(t)}{dt} \neq f(t)$$

Except in the specific case where $f(t)$ is a constant.

The correct way to generate a sine wave with time-varying frequency $f(t)$ is thus with:

$$y(t) = \sin (2 \pi \int_{0}^{t} f(\tau) d\tau)$$

In synthesizer-speak the register/variable accumulating the instantaneous frequency to evaluate the $\int_{0}^{t} f(\tau) d\tau$ quantity is called a phase accumulator.

http://en.wikipedia.org/wiki/Frequency_modulation

http://en.wikipedia.org/wiki/Instantaneous_phase


How can you perform echo cancellation on features of audio signals? I am working on a project that takes some features of audio signals and classify them with a neural network. However, the audio signals are added with their own echo with the environment. What would be some of the processing techniques that can be applied on the features, including chroma stft and Mfcc.

well, while applying a room impulse response to audio (which is what echoes are) is a linear operation on the time-domain signal, it's generally not on features.
So, you're looking to find a non-linear mapping of a somewhat high-dimensional feature vector space to a "cleaned up" feature vector space with some well-known metric on how well that mapping performs as cleanup....
That sounds exactly like the kind of job that neural networks are good at. And that's not a coincidence: When you look under the mathematical hood of all your ML algorithms, there's the Universal Approximation Theorem that says that a NN that's sufficiently large can represent any non-linear mapping.
So, really, that's a hard problem to do analytically, and it's thoroughly reasonable to assume that "echo removal on MFCC features" can most sensibly be done with a Neural network, most likely even the same that you directly train to classify your audio signals.


what is the maximum instantaneous frequency of Dirac delta? I wish to find the maximum frequency of $\sin(\varphi(t))$, where $\varphi(t)=u(t)$. We have $\omega(t)=d\varphi(t)/dt=\delta(t)$. What is the the maximum instantaneous frequency $\omega_{\max}$?

It's not very sensible to even try to assign an instantaneous frequency to signals like yours.

Your specific example fails for more fundamental reasons than usual however, namely for the fact that the Dirac delta is not a function but something mathematicians call a "distribution".

This for example implies that we cannot assign a value to $\delta(0)$ even though we know $\delta(x)=0$ for all $x\neq0$. The definition of the delta distribution is, that $\int_\mathbb{R} \delta(x) f(x) dx=f(0)$. From here it follows that $\delta(0)$ cannot have a finite value. However, stating that it is "infinite" there is also not enough. There's simply no way to think of any value at $0$ that makes sense.

So the answer to your question is that the instantaneous frequency is $0$ almost everywhere, and undefined at $t=0$. So the maximum (and minimum!) value of $\omega$ is $0$!

You can avoid this situation by replacing the unit step with a smooth approximation of a unit step. With this you can give your maximum instantaneous frequency any value you desire, even if you take a limit of a function sequence that results in the unit step itself. So the answer is absolutely arbitrary.


Does Function to Function Multiplication and Convolution, Follow Laws of Convolution Apologies this is long – I do not want to leave any ambiguity.

I have a hopefully simple question, that does not seem to be covered by the DSP books I have in regards to theory of operation.

I will use the ‘*’ as the operation for convolution, and ‘x’ to denote multiplication. For multiplication, the operation is an element by element multiplication, as opposed to a constant multiplication.

Assume that there is a continuous sine wave, which will be sampled at 192kHz, with frequency 1kHz and amplitude 1, which lasts for 1 second, denoted as sin(n).

We also have the rectangular function, with the same sample rate, which lasts for 5ms, about the centre of the 1 second interval, denoted as rect(n).

There is a filter function that based on the sample rate has a cut off at 24kHz, is defined as symmetrical, linear phase, and has 501 taps – denoted as h(n).

What I want too know is, are the two operations described below valid in terms of LTI systems, convolution, and multiplication, and the distributive/associative/commutative laws.

Process A). The sine wave is multiplied by the rectangular function to produce a truncated sine wave, which starts and ends at the zero crossing of the sine wave. Essentially creates a gated sine wave. The resultant truncated sine wave is then applied to a filter – convolution.

 rect(n) x sin(n) = A(n)    →     A(n) * h(n) = Q(n)

Process B). The sine wave is convolved with the filter, and the rectangular function is convolved with the filter, and the result of the two outputs, are then multiplied together.

sin(n) * h(n) = B(n),       rect(n) * h(n) = C(n),    B(n) x C(n) = R(n)

The question is, is Q(n) = R(n) ?.

I used Octave to implement this, and they are not the same, but the errors are small, and visually the Process A output looks different to Process B output about the rectangular function discontinuities.

[Reason for asking this question : it has been claimed -

You have two input functions:
1) a sine
2) a boxcar function
The product of the two gives you the gated sine. You present this to a steep linear phase filter. This filter then rings, not because of 1), but because of 2). 1) occupies a single spectral line,but 2) occupies an infinite spectrum.
As the filter is a low-pass and linear phase it spreads function 2) to the left and to the right of the time interval it originally occupied. The sine itself is not affected by the filter because it resides entirely in the passband. So if you then cut out the sine part and join the two ringing parts then you are looking at the ringing as caused by a boxcar.]

Thanks and Regards,
Code_X

In a nutshell, you are asking whether the following statement is true:

$$[x(n)\cdot y(n)] *h(n) = [x(n)*h(n)] \cdot [y(n)*h(n)]$$

In words, if the convolution operator is distributive with multiplication.

The short answer is no.

An easy example would be to think of two sequences $h(n)$ and $x(n)$ such that their convolution is $0 \ \forall n$ (if you can't think of any, check this out). We can see that in this case the right side of the equality given above would be $0$ no matter what $y(n)$ is, but the left one would not (at least not for all $y(n)$).


Problems Using FFT to Compute Impedance in a Model Neuron I'm a neuroscientist currently investigating the resonance properties of a single neuron model that a colleague and I have constructed. The language we code in is Julia, which I hope is similar enough to Python for most people to understand. I've run in to a very particular problem, so bear with me please.

In order to investigate the resonance properties of the model neuron, I am applying an exponentially swept sinusoid current to it and analyzing its voltage response. In my field, such a current is called a ZAP current, though a more familiar term is probably Chirp signal. Just so there is no confusion, the current is generated from the following function \begin{gather*} ZAP(t) = I_{app}sin(2\pi f(t)t)\\f(t) = (0.001)\frac{f_{min}}{Lt}(e^{Lt}-1)\\ L = \frac{\ln{\frac{f_{max}}{f_{min}}}}{t_{max}} \end{gather*} where $t_{max}$ is the time duration of the simulation, $f_{min}$ and $f_{max}$ are the starting and ending frequencies of the sweep, and $I_{app}$ is the magnitude of the sinusoid. The voltage response of the model neuron to this input looks like this:

Clearly, there is a resonance in the model at whatever the input frequency is at around 50,000 ms. This means there should be a peak in the impedance of the model at this frequency. The impedance is calculated as the ratio of the discrete Fourier transform of the output voltage to the input current $$ Z = \frac{FFT(V(t))}{FFT(ZAP(t))} $$ I then plot the magnitude of the impedance against the natural logarithm of the frequency. I expect to get something like the one on the left, but instead I obtain the one on the right:

As I vary the parameter $E_K$, the resonance frequency changes, and I expect to see this reflected in the impedance-frequency plot, but every one of the impedance-frequency plots looks almost identical to the one shown above despite the voltage responses looking quite different.

I have spent hours going through my code line by line and experimenting with different modifications to figure out where I'm going wrong, but have not been able to produce the desired output. Here is the stripped down code:

## Time Parameters
dt      = 0.01         # ms
tmax    = 100000        # ms
tspan   = convert(Int64,tmax/dt)
t       = linspace(dt,tmax,tspan)

## Simulation and Analysis
# magnitude of sinusoid
Iapp = 1.0
# minimum frequency
f0 = 0.1
# maximum frequency
ff = 15.0
# create an array of the frequency sweep values for plotting
Zfreq = ZAPfreq(f0, ff)
# optimize Fast Fourier Transform for an array of size and type ZAP()
P = plan_fft(ZAP(Iapp, f0, ff))
# calculate FFT for input ZAP()
inFFT = P*ZAP(Iapp, f0, ff)
output    = simulate(ZAP(Iapp, f0, ff))
# calculate the FFT of the membrane voltage response
outFFT = P*output
# calculate the impedance of the cell at each frequency ZAPfreq[i] from Z = V/I
Z = outFFT ./ inFFT
# calculate the magnitude of the complex impedance
MagZ = abs(Z)
# center the DC response (f = 0 at center of frequency axis)
shiftZ = fftshift(MagZ)
# plot the impedance at each frequency
plot(Zfreq, shiftZ, grid = false,
  legend = false,
  linecolor= :black,
  xscale= :ln,
  xlabel = "log(\$f\$)",
  ylabel = "Impedance (Ω)",
  title = "Kispersky-Caplan Impedance Response with \$E_K\$ = $EK",
  top_margin=5*mm,
  bottom_margin=5*mm,
  left_margin=5*mm,
  right_margin=10*mm,)

Here, the $\texttt{simulate()}$ function is the one that actually carries out the simulation of the model. Any help in determining the problem is greatly appreciated.

I don't know how much background you have with FFT, so I apologize if I list a few things you already know.

In general, I think people have a tendency to throw an FFT at a problem and hope it works, but you may need to consider what the transform is really all about. First, when you take an FFT, it's implicit in the calculation that your signal is periodic. So, even though you're feeding it a finite number of samples with frequency content only between f0 and ff, the calculation makes the assumption that before time zero and after your samples, the waveform repeats forever. If your waveform's beginning and end don't match perfectly, there's going to be a lot of frequencies you wouldn't expect in your FFT. Basically, it's like there's a discontinuity. So, often, you'll window the signal to force it to zero at both ends so each periodic clone of your waveform meets the next smoothly. This obviously also has consequences for the frequency content of the FFT, so you need to do this carefully also.

Second, if your waveforms have any DC content (a constant offset, frequency = 0 Hz), it usually makes sense to remove it before your FFT unless DC is of interest to you. Often, low frequencies are huge and swamp the stuff you care about. Doing some high-pass filtering with a very low cutoff will help you see the interesting content of your plot more easily.

Third, when you view your plot on a linear frequency axis, you'll see that there are Fourier coefficients for both positive and negative frequencies, and the plot is symmetrical about 0 Hz. For qualitative assessments, you probably only need to look at the positive half.

Fourth, consider that you're taking the ratio of output over input, and anywhere your input FFT is near zero, that's going to blow up, and may make it hard to see what's going on.

So, I'd say, do some preprocessing of your signal to remove things you don't care about, and to prevent edge effects. Then post another plot and let's see how it looks. Since you're simulating the input and output here, there may not be a need, but in practice, I'll often use Welch's method to compute power spectral densities, and then compute an H1 or H2 estimation of the frequency response. Here, we can hopefully get a result with a simple analysis.


How to remove noise using mean and median filters? I can remove the noise using moving average filter. The MATLAB code is shown bellow.

But how can I remove the noise using mean and meadian filters?

You do not need loops. This implements a moving average filter, a median filter, and a median filter followed by a moving average:

%x = noise1;
x = randn(1,2000)+1;
lengthMean = 5;
lengthMedian = 7;
xFiltMedian = medfilt1(x,lengthMedian);
xFiltMean = filter(ones(lengthMean,1)/lengthMean,1,x);
xFiltMedianMean = filter(ones(lengthMean,1)/lengthMean,1,xFiltMedian);
figure(1);clf;hold on
plot(x(500:1500),'b');
plot(xFiltMedian(500:1500),'r');
plot(xFiltMean(500:1500),'g');
plot(xFiltMedianMean(500:1500),'c');

As said by @jojek the median filter is often more efficient than an averaging filter for spikes and outliers. And if you choose to apply both filters, it is generally more useful to apply the  median first (as done in the code above).

The choice of filter length should normally be related to properties of your signal and the noise (spectrum, distribution), and to the goal  you are assigned with. But the length of the signal is usually not a criterion.  For symmetry  properties, odd-length filters are practical.

Without further details, you can start with longer filters (e.g. lengthMedian= 37), see how they perform (how they reduce the  noise but also degrade the signal's properties), and reduce the length progressively to the equilibrium point, where shorter filters do not yield better results.

However, the noise in the  valley suggests it is not stationary, so you may need more adaptive filters (e.g. with a varying length), or a little more involved like a weigthed median filter that limits the location and flattening issues you get with long median filters. It is discussed on SE Cross validated.


Recommendation for book - Writing DSP code in C I am looking for some good book, that simply show how you actually write a code in C, to do all the main DSP methods .

FFT.
Low-pass and high-pass filters.
Auto-correlation.
Noise processing.

And all the basics of DSP, from theory into a real code in C.

For example, I have got 1000 samples, now i would like to compute FFT of it, remove noises, then bring back to the time axis.

Is there something good that cover all of these?

I can recommend you two books about DSP for C language.

  Embree P. M. - C Language Algorithms for Digital Signal Processing

It is old and you can easily get it second-hand for a decent price. It covers pretty much all 4 topics that you described.

The other one I recommend is:

  Malepati H. - Digital Media Processing: DSP Algorithms Using C

It covers way more topics (including error correction algorithms and image/video processing algorithms). The nice thing about that one is that test signals are provided to test your implementation.

I just recalled that there is another book on my desk which has some C code for DSP in it.

  Press W. H. - Numerical Recipes: The Art of Scientific Computing

There are two chapters about the Fourier Transform and its applications.

One last suggestion - if you can, please use an existing C library instead of writing your own tools. There is no point in re-inventing the wheel.


How to interpret a spectrum with FFT mirroring? I have this code which interprets a few data points on magnetic field strength into a spectrum, as shown below:
c = fft(magfield);
c = c*2/n;
ca = abs(c);
dt = 1/fs;
f = (0:n-1)/(500*dt);

for i = 1:length(ca)
    if ca(i) > 0.002608
        ca(i) = 0;
    end
end

spectrum = fftshift(ca);
plot(spectrum)

Where fs is 5e5 and n is 50000000. what i get is a graph which looks like the one below:

I have no idea how to interpret the spectrum, where the peak frequencies are and why this mirroring effect is happening. Could you guys please shed some light on what's going on?
Thanks!

Your field data is real (as in "not complex") and hence your spectrum is conjugate symmetric.
fftshift() rotates the spectrum so that 0 Hz is in the middle of the vector. Your x-axis on the plot should be fs*[-nFFT/2:nFFT/2-1]/nFFT; where $nFFT$ is the FFT length and $fs$ is the sampling rate.
Remember that time discrete signals are periodic in frequency the second half of the spectrum represents the negative frequencies.


How to transform or distort or skew a rectangular image to fit arbitrary quad I have a rectangular image which need to transform so the image fits an arbitrary convex quad. Easy? Except I have an important constraint, the features on the original image must maintain their relative area after the transformation.

I've tried a perspective transformation, but of course feature's relative area's change significantly with perspective.

The (hand-drawn) image below illustrates the effect I need. The original is on the left and the transformed is on the right. Notice that all 4 features maintain same relative area.

I would be happy with solution that maintained the areas approximately.

Assume that the center of the square maps to the center of the vertexes.
The vertices in the output are known. Assume they are $(x_1,y_1),...,(x_4,y_4) $
Their center is $ ( \frac{x_1+x_2+x_3+x_4}{4},\frac{y_1+y_2+y_3+y_4}{4}) $
Now you have 4 degrees of freedom of where to put the points where the color change.
$ ( w_1 x_1 +(1-w_1) x_2 ,w_1 y_1 +(1-w_1 ) y_2 )  $, and so on

Areas of each of the 4 parts equal to the areas of two triangles (Can be calculated analytically)
All of these together should give you enough constraints to find $w_1,w_2,w_3,w_4$
Once you found the weights, you can create a grid based mapping that transforms squares into the required forms.


Why in FSK the phase shift is allowed? FSK modulation encodes information into the frequency of the signals. But when I checked the reference books, the formal expression for a frequency modulated signal is:
$x_i(t) = \alpha cos(2\pi f_it + \theta_i)$
Which means, for a binary FSK, we can consider $x_1(t)$ and $x_2(t)$ with different frequencies $f_1$ and $f_2$ (which make sense) and different phase $\theta_1$ and $\theta_2$ (which is odd to me), respectively.
(1) My first question is, why should we allow the signals have different phases? I cant see a reason for this and as far as I know, different phases result in phase discontinuity of the transmit signal.
(2) And the second question is, how this difference in phase is related to coherent and non-coherent receiver?

My first question is, why should we allow the signals have different phases?

Counterquestion: Why should they not be allowed? Enforcing continuous phase can make systems more complex, without any benefit.

I cant see a reason for this and as far as I know, different phases result in phase discontinuity of the transmit signal.

Depending on how you do that transition, that's not necessarily the case: you can pulse-shape $\theta_i$, or you can filter $x$.
It's less often a problem than old textbook make it out to be, to be honest.

(2) And the second question is, how this difference in phase is related to coherent and non-coherent receiver?

Not at all. What the receiver does says nothing about what the receiver does!