system bit rate vs channel coding We are  usually asked about the full system bit rate. However, we use channel encoding in practical systems, for example conventional coding with rate $(1/3)$. My question is about the full system bit rate in such system which employs channel coding. Should it be calculated including the coding bits or after extracting the code bits? I mean if we have system whose bit rate of $12bps$ and employs channel encoding of rate $(1/3)$. Should the bit rate be $12bps$ or $4bps$ since the other $8$ bits are used for coding?

Both indicators are used. It depends what you mean in speaking of your indicators.
One of common usages is to compare systems where we use indicators of the same type. For example,

to compare 5G NR and LTE RANs, we usually use max throughput which is defined as the multiplication of number of antennas, number of subcarriers and QAM order per subframe. There is no error probability here.
To compare channel decoders e.g. iterative belief propagation LDPC decoder and Ordered Statistic decoder, it does not make sense anymore if error probability is ignored and, therefore, the code dimension (number of bits before encoding) must be used.


Determine impulse response given input and output: which ROC? Let's suppose I have to find the impulse response of a discrete time LTI system given a specified input and its output through the system. I think I'm going to get the $\mathcal Z$-transform of input-output and then as everyone knows
$$H(z)=\frac{Y(z)}{X(z)}$$
The problem is: I want it in time domain, so I need to inverse transform $H(z)$. How will I choose the ROC (hence causality / noncausality)?

The function $H(z)$ without a specified ROC is generally not sufficient for uniquely determining the time domain sequence $h[n]$. You need more information about the desired $h[n]$. There are three types of ROCs:

the ROC is an annulus ($R_1<|z|<R_2$): $h[n]$ is a two-sided sequence (non-causal)
the ROC is defined by $|z|>R$: $h[n]$ is a right-sided sequence
the ROC is defined by $|z|<R$: $h[n]$ is a left-sided sequence

If you're looking for a stable system, then the ROC needs to include the unit circle $|z|=1$. If that does not uniquely define the ROC, then you need to know if you're looking for a left-sided, a right-sided, or a two-sided sequence.


Problem getting IFFT to work I've run into an issue fiddling with the FFT/IFFT routines given here (http://www.dspguide.com/ch12/3.htm). I've rewritten the code to C++ and played with it. The idea is to read a raw wave file (encoded as 32 bit floats),
find its FFT (complex) and then use that FFT as input to the InverseFFT function.

I expected to get the original sound data (samples) and I do for FFT lengths of up to 2048. For greater lengths (>= 4096 ), the code below doesn't (seem to) work. As you can see, my test consist solely of checking values
at indexes 20,30 and 40 for both the original samples and the result (inverse FFT) array. I've tried this on several files (created with Audacity) with the same results.

    int main(int argc, char *argv[])
{
    unsigned fs;
    float *Ix;
    QCoreApplication a(argc, argv);
    cout << "opening file...";

    char *pData= OpenRaw("C:\\test\\250_350_c.raw",fs);

    float* Rx = (float*)pData;
     cout << "original data { [20,30,40] } " << Rx[20] << " " << Rx[30] << " " << Rx[40] << endl;
    const unsigned numSam = 2048;
    Ix=new float[numSam];
    for(int i=0;i<numSam;i++)
        Ix[i]= 0;

    FFT(Rx,Ix,numSam);

    float *Ifft = InverseFFT(Rx,Ix,numSam);

    cout << "IFFT data { [20,30,40] } " << Ifft[20] << " " << Ifft[30] << " " << Ifft[40];
    return a.exec();
}
char* OpenRaw(const char *fileName, unsigned& fs) {
    char *pData;
    unsigned fileSize = 0;
    int bRead = 0;

    ifstream file (fileName, ios::in | ios::binary);
    if (file.is_open())
    {
        cout << "reading file..." << endl;

      file.seekg (0, ios::end);
      fileSize = file.tellg();
      cout << "file size " << fileSize << endl;

      pData = new char[fileSize];
      file.seekg (0, ios::beg);
      bRead = file.readsome(pData, fileSize);

      if(file.fail())
          cout << " error reading file..." << endl;
      else
          cout <<  "bytes read.." << bRead << endl;
      file.close();

    }else{
        cout << "error reading file";
    }
    fs = fileSize / sizeof(float);
    return pData;

}
float *InverseFFT(float *Rx, float *Ix, int fftSize){
    float *pWaveData = new float[fftSize];
      for(int i=0; i<fftSize; i++){
          Ix[i] = -Ix[i];
      }
      FFT(Rx,Ix,fftSize);//calling the FFT function
        for(int i=0; i<fftSize;i++){
         pWaveData[i] = Rx[i]/fftSize;
         Ix[i] = - Ix[i]/fftSize;
        }
    return pWaveData;
}

I haven't posted the FFT code here because I've used that code in the past and it always worked OK for me (I can still post it here if someone wishes to experiment with the code).

Does anyone know what might be causing this behaviour?

UPDATE 1: I've implemented the second method found in the link given by Peter below, but the result is still the same (the IFFT works for fft lengths of only up to 2048). My original IFFT method (given above) is actually the method #4 found in that same link.

float *InverseFFT(float *Rx, float *Ix, int fftSize){

FFT(Rx,Ix,fftSize);
float *pWaveData = new float[fftSize];

  for(int i=1; i<fftSize/2; i++){
      pWaveData[i] = Rx[fftSize-i]/fftSize;
      pWaveData[fftSize-i] = Rx[i]/fftSize;
  }

return pWaveData;

}

UPDATE 2: After googling some more, I found this ( https://www.mathworks.com/matlabcentral/newsreader/view_thread/247760?requestedDomain=www.mathworks.com ).
Well, unfortunately, it seems that longer fft lengths may adversely influence an FFT/IFFT algorithm, especially if it's a power of 2 based algoritm.

Thank you

I solved this by using a different FFT implementation, one that is using a mixed radix algorithm under the hood ( the project is avalaible here https://github.com/itdaniher/kissfft ). The KissFFT supports both a complex as well as a real FFT so I'll post code for both. The real version is up to 45% faster than the complex one. Note that in order to get the original samples, the results of the IFFT have to be scaled by fft size (fft size doesn't have to be a power of 2 but for the real FFT it must be an even number).

The complex FFT/IFFT

    #include "kiss_fft.h"
    #include "kiss_fftr.h"

        int main(int argc, char *argv[])
            {

                unsigned fs;

                QCoreApplication a(argc, argv);
                cout << "opening ...";
                char *pData= OpenRawA("C:\\test\\250_350_c.raw",fs);

                const unsigned numSamples = 4096;//this doesn't have to be a power of 2
                float* Rx = (float*)pData;

                       kiss_fft_cpx *fIn = new kiss_fft_cpx[numSamples];
                       kiss_fft_cpx *fOut = new kiss_fft_cpx[numSamples] ;
                       kiss_fft_cpx *fIfft = new kiss_fft_cpx[numSamples];
                       //populate the (complex) input buffer
                        for(int i=0; i<numSamples;i++){
                            fIn[i].r = Rx[i];
                            fIn[i].i = 0;
                        }

                        cout << "original data { [20,30,40] } " << fIn[20].r << " " << fIn[30].r << " " << fIn[40].r <<  fIn[50].r << " " << fIn[60].r <<endl;
                        kiss_fft_cfg cfg = kiss_fft_alloc( numSamples ,0 ,0,0 );//configure the FFT function to use a forward FFT

                        kiss_fft( cfg , fIn , fOut );//calculate FFT

                        cfg = kiss_fft_alloc( numSamples ,1 ,0,0 );//configure the FFT function to use an inverse FFT

                        kiss_fft( cfg , fOut , fIfft );//calculate inverse FFT

                       free(cfg);

                cout << "IFFT data { [20,30,40] } " << fIfft[20].r/numSamples << " " << fIfft[30].r/numSamples << " " << fIfft[40].r/numSamples << " " << fIfft[50].r/numSamples << " " << fIfft[60].r/numSamples;
//delete allocated resources
                return a.exec();
            }

The real FFT/IFFT

 int main(int argc, char *argv[])
    {

        unsigned fs;

        QCoreApplication a(argc, argv);
        cout << "opening ...";
         char *pData= OpenRawA("C:\\test\\250_350_c.raw",fs);

        const unsigned numSamples = 4096;
        float* Rx = (float*)pData;

               kiss_fft_scalar *fIn = new kiss_fft_scalar[numSamples];
               kiss_fft_cpx *fOut = new kiss_fft_cpx[numSamples] ;
               kiss_fft_scalar *fIfft = new kiss_fft_scalar[numSamples];
               //populate the input buffer
                for(int i=0; i<numSamples;i++){
                    fIn[i] = Rx[i];

                }
                cout << "original data { [20,30,40] } " << fIn[20] << " " << fIn[30] << " " << fIn[40] <<  fIn[50] << " " << fIn[60] <<endl;
                kiss_fftr_cfg cfg = kiss_fftr_alloc( numSamples ,0 ,0,0 );//configure the FFT function to use a forward FFT

                kiss_fftr( cfg , fIn , fOut );//calculate FFT

                cfg = kiss_fftr_alloc( numSamples ,1 ,0,0 );//configure the FFT function to use an inverse FFT

                kiss_fftri( cfg , fOut , fIfft );//calculate inverse FFT

               free(cfg);

        cout << "IFFT data { [20,30,40] } " << fIfft[20]/numSamples << " " << fIfft[30]/numSamples << " " << fIfft[40]/numSamples << " " << fIfft[50]/numSamples << " " << fIfft[60]/numSamples;
//delete allocated resources
        return a.exec();
    }


Sampling Theorem and Signals Explained to a Mathematician Let $f:\left(-\frac T2,\frac T2\right)\to\mathbb{R}$ for some $T>0$. The Fourier coefficients of $f$ are $$\left\{\begin{matrix}a_0&=&\displaystyle\frac 1T\int_{-T/2}^{T/2}g(t)\;dt\\a_k&=&\displaystyle\frac 2T\int_{-T/2}^{T/2}g(t)\cos\left(k\frac{2\pi}Tt\right)\;dt\\b_k&=&\displaystyle\frac 2T\int_{-T/2}^{T/2}g(t)\sin\left(k\frac{2\pi}Tt\right)\;dt\end{matrix}\right.\;\;\;\;\;(k\in\mathbb{N})$$
and the Fourier polynom of degree $n\in\mathbb{N}$ itself is $$\mathcal{F}_n[f](x):=a_0+\sum_{k=1}^n\left(a_k\cos\left(k\frac{2\pi}Tx\right)+b_k\sin\left(k\frac{2\pi}Tx\right)\right)\tag{1}$$ Now, for any $g\in L^2$, where $$L^2:=\left\{g:(-\pi,\pi)\to\mathbb{C}:g\text{ is continuous almost everywhere}\right\},$$ the complex Fourier coefficients $c_k$ of $g$ satisfy Parseval's identity, i.e. $$\sum_{k\in\mathbb{Z}}\left|c_k\right|^2=\frac{1}{2\pi}\left\|g\right\|_{L^2}^2\tag{2}$$

Now, we can interpret $f$ as a $T$-periodic signal that we want to transmit over a channel. I'm asking myself the following questions:

What is the signal bandwidth?$^1$
What is the frequency spectrum?
What is the channel bandwidth?
What does the sampling theorem state?
I've read that the square of the $k$-th coefficient, i.e. $a_k^2+b_k^2$, is proportional to the "energy contained in this harmonic". I think this is somehow related to $(2)$. But I have no idea what is meant by energy.

I've read all the definitions, but (as I'm a computer science and mathematics student with almost no relation to physics) I wasn't able to answer the questions above by myself.

I know that $(4.)$ has something to do with the fact, that we need to use some kind of discrete Fourier transformation, i.e. we consider an equidistant grid $$x_j=-\pi+j\frac \pi N\;\;\;\text{for }j=0,\ldots,2N$$ and (maybe) use the composite trapezoidal rule to approximate $$c_k=\frac{1}{2\pi}\int_{-\pi}^\pi f(x)e^{-ikx}\;dx\approx\frac{1}{2N}\sum_{j=0}^{2N-1}f\left(x_j\right)\exp\left(-ik\frac{2\pi}Tx_j\right)$$ So, I think the sampling theorem targets the question how huge do I need to choose $N$ in order to reconstruct (what exactly does that mean? If the Fourier polynom has infinitely many summands we cannot reconstruct the original function) the signal and at which size of $N$ it would be pointless to choose an even bigger $N$.

$^1\;\;$My lecture notes distinguish whether or not $\displaystyle\mathcal{F}[f]:=\lim_{n\to\infty}\mathcal{F}_n[f]$ is actually finite (i.e. $\exists k_0\in\mathbb{N}:\forall k\ge k_0:a_k=b_k=0$). They state, that "the signal bandwith is the difference between the lowest and the highest frequency"?  What does that mean? In the first place, I wasn't even sure if signal bandwidth is a property of $f$ or $\mathcal{F}[f]$. Considering $f$, I would not understand why $f$ should have multiple frequencies and how I would calculate them. So, I think we need to consider $\mathcal{F}[f]$. More precisely, I think we need to look at the "frequencies" of the sine and cosine functions in the $k$-th summand in $(1)$. Obviously, these frequencies are both equal to $$\mathcal{f}_k:=\frac kT$$ and independent of $f$. But, I still don't understand how I can determine the "highest" frequency.

"Highest frequency" is shorthand for the highest k/T for which a and b sub k are non-zero.

Considering f, one reason it might consist of the sum of multiple sinusoidal components is that this is a solution to many ordinary linear low-order differential equations.


How to calculate PSD if i have power vs freq spectrum? i have a signal something like this

how to calculate PSD?
i have freq from 191000200 to 196500000
time gap or delta t = 1 or 2 sec, (what i mean is [tn] - [tn-1])
i will be computing in python any suggestion from scipy or any lib will be helpful
thanks in advance

i used this formula saw on wiki, https://en.wikipedia.org/wiki/Spectral_density#Power_spectral_density
it works okay for me,
thanks for the comments, it helped me to understand the data what I was receiving.


What can we deduce about variance when we are given the noise spectral density? Given the noise spectral density ${N_0\over 2} = a$, what can be deduced about the variance?

In the case of a two-symbol canal(/cable), and when $T_0 = 1$.

Is it that, $\sigma_X^2 = {N_0 \over 2} *T_0$?

Given a low-pass signal $x$ that is a realization of random variable $X$ with constant spectral density $a=N_0/2$ in $f=-B\ldots B$ (and zero otherwise) we can calculate the mean power of that signal by $\bar P = 2aB=N_0B$, where $B$ is the signal bandwidth. The mean power of a random variable is defined as the expectation of $X^2$, i.e. $\bar P = \mathrm E[X^2]$. The correspondance between mean power and variance $\sigma_X^2$ of a random variable $X$ is
$$
\sigma_X^2=\mathrm E[X^2]-m_X^2
$$
where $m_X$ is the mean value of $X$. That is, for zero mean (which is often true for noise) we get
$$
\sigma_X^2=\mathrm E[X^2] = N_0B
$$

Concerning your second question, I assume that $T_0$ is the symbol duration. Thus the symbol rate $R_\mathrm{S} = 1/T_0$. The bandwidth of a signal with symbol rate $R_\mathrm{S}$ is often approximated by $B=R_\mathrm{S}/2$ and thus
$$
\sigma_X^2 = N_0B = N_0R_\mathrm{S}/2 = N_0/(2T_0)
$$
(You can also check the units in your equation: $N_0$ is power per Hz and $T_0$ is in sec. But the result should be a power)


Avoid Aliasing when filming a moving car I want to film a moving car with a rolling circumference: 190,5 cm and wheel circumference: 60 cm and 5 spokes.

I record with 24Hz. What can be the maximal speed of the car to avoid aliasing?

Let's say $\omega_w$ is wheel rotational speed, $v$ is car speed, $r = 0.3m$ ($2\pi \cdot 0.3 = 1.88m = $ rolling circumference?) is the radius of the wheel and $f_s = 24Hz$ is the sampling frequency.

Then you need to sample an event with a period of $5\omega_w$ (stokes are indistinguishable from one another). Thus, $f_{max} = \frac{5\omega}{2\pi}$.

Then, by Nyquist:

\begin{equation}
f_s \gt 2f_{max} = \frac{2 \cdot 5\omega_w}{2\pi} = 5 \cdot \frac{v}{r \cdot \pi}
\end{equation}

Therefore:

\begin{equation}
v \lt \frac{fs \cdot r \cdot \pi}{5} = \frac{24 \cdot 0.3 \cdot \pi}{5} \approx 4.52 m/s \approx 16.28 km/h
\end{equation}


Finding Homography Matrix Using Lines I have some Points and Lines from original image and the corresponding points and lines from top view image.

-> I tried to find homography using points

Used OpenCV function findHomography(objectPointsPlanar, imagePoints) :

It takes 4 points from original image and 4 points from top view image(required image) and it returns homography matrix of 3x3.

Used mathematical formula to find homography matrix (Didn't use any library) :

Please look at the attached images for the formula to find homography matrix using 4 points by your own.

I have tested the formula by writing it in c++. It gives the homography matrix same as given by opencv function findHomography().

PS : I used the determined homography matrix to warp the image, but I am not happy with the output. It doesn't come to be accurate.
So I am thinking to use lines to find homography matrix.

-> Homography determination using Lines :

I tried to look for any library which can find homography matrix using lines, but I didn't get.
I also looked at research papers to get the mathematical formula to find the homography matrix, but was not able to understand the formula.
I am sharing the link of a research paper  : A new normalized method on line-based homography estimation

So could someone please help me to find the homography matrix using lines.

Thanks in advance,
Regards.

I think you can't do Homography on lines.

I don't remember the equations but I don't think lines in one image are kept as lines in the other image.

What you can do is treat a line as set of points and have Homography estimation using many point samples.


Nyquist frequency - Sampling theory If a signal contains frequencies from $500$ to $1000\textrm{ Hz}$, what is the Nyquist sampling frequency, $2$ times $500$ or $2$ times $1000$?

Twice the highest frequency of the signal, or twice the bandwidth of the signal?

Both.  For baseband sampling (thus using only low-pass filtering), a sampling rate over over 2x the highest frequency is necessary.  For band-pass under-sampling, sampling at a bit over twice the bandwidth will work, as long as you remember the original band and don't fold (alias) the given band onto itself.  1000 Hz sampling will fold and invert a (slightly narrower than) 500 to 1000 Hz signal down to 500 to 0 Hz, and thus the signal won't be aliased with itself, but with an empty frequency band (after band-pass filtering stuff out, if needed).  Thus, no information will be lost (aliased), as long as you remember the original frequency band (e.g. it was not baseband).  Reconstruction may require a band-pass anti-aliasing filter, rather than a low-pass anti-aliasing filter.  Using knowledge of the original band, you can de-invert and upconvert the samples back to the original lower side-band of a 1000 Hz carrier, for example.


How to find a know signal in another signal which contains other signals along with the signal i need in python I am working on BPSK with sound waves a carrier waves.

Actually I am modulating my data with some known sync bits along with the information I need to send. So the things I will be having at demodulation side is these sync bits and carrier wave.

So I am modulating these sync bits with carrier and cross correlating this modulated (over sync bits) wave with received wave at the transmiter. By doing this I want to find the exact time position (where these sync bits appear in received wave).

For this I am first crosscorelating received wave and my sync wave then, when I plot it I wil surely get a maximum where the sync wave in both signals coincide.

But I am not sure if this maximum is the position I need (when there is a DC offset or large amplitudes in the received signal, it can also create an extreme value)

So to overcome this, I autocorelated my sync wave with itself and took the
absolute difference between the all the values of crosscorelated matrix and
autocorelated value. Then I assumed which ever value has the least value as the
synced location.

This is also not working (I think as iam not sure that
the location iam getting is midpoint of sync or starting of sync). Is this
approach correct?

Or are there any other approaches to find the starting sync position in received wave?

I also needed help in knowing whether if we use the first method (wher evr maxima comes in crosscorelation taking tht as sync location) is this location I got, is sync starting location or midpoint of sync?

I made a little research on detecting the sync_start location on my received signal.

First I modulated the bits (sync_start bits + data bits + sync_end bits) with a sine signal as my carrier, with a frequency of $20 \text{ kHz}$ and a total length of 10000 bits (including both sync_start and sync_end).

I then converted this modulated wave into a .wav file
and played it on my PC, while recording simultaniously on an another PC, using Pyaudio (Python). That gave me a recorded .wav - my received signal.

Now as I knew my sync bits value, I modulated them on the same carrier wave and cross correlated this with the received wave.

Case 1:

$\text{bit length} = 10000$ (including sync_start and sync_end),

$\text{number of sync_start bits} = \text{sync_end} = 2287$

The plot of cross correlation is as shown :

Case 2:

$\text{bit length} = 15000$

case 3:

$\text{bit length} = 20000$

As the data length increases(keeping our sync_start and sync_end same in all cases), the exact time stamps of sync_start_mid index and sync_end_mid index gets more accurate.

Note : For data bits here, I took all of them as $1$'s -- that might be the reason for the smaller amplitudes in the middle. I will add random data and see if I am getting the same inference.

Please correct me if there is any mistake in my inference. As far as now I am now able to exactly pinpoint the synchronization between transmitted and received waves.

I also wanted to know the limitations of length of the sync bits (because when my sync bits length gets smaller, the dominant peak location cannot be found from cross correlation).

And I am using a BPSK scheme.


Should the phase coefficients of the DFT of a real, even input signal all be zero? The phase coefficients of a real, even input signals should all apparently be $0$ or a multiple of $\pi$. That's a property of the DFT I've learned about in the Audio Signal processing course on coursera. However, the example they use completely confuses me.

They create a triangle wave that looks like this:

and then showed the phase spectrum to be this:

They explained that this phase spectrum was not zero, because the original triangle function wasn't even because it wasn't centered around zero - it had a phase shift that was messing things up. So they adjusted the positioning of the triangle wave so the peak was at zero, like this:

And the phase spectrum duly went to zero (note the minuscule scale on the $y$-axis)

Bearing in mind I'm new to DSP, the thing that is really confusing is that surely both forms of the triangle wave are even - shifting by half a wavelength just creates another even function, doesn't it? Why must the triangle wave be positioned so the peak is at the center rather than the trough, or vice-versa?

You are right, of course it's not necessary to position the peak of the triangle in the center. The only condition for the DFT of a sequence to be real-valued is that the sequence satisfies

$$x[n]=x^*[N-n],\qquad 0\le n<N\tag{1}$$

where $*$ denotes complex conjugation, and $N$ is the DFT length. If $x[n]$ is real-valued, $(1)$ reduces to

$$x[n]=x[N-n],\qquad 0\le n<N\tag{2}$$

Conditions $(1)$ and $(2)$ require the periodic continuations of $x[n]$ (period $N$) to be even with respect to $n=0$.

From your first figure it appears that your signal satisfies

$$x[n]=x[N-1-n]\tag{3}$$

with $N=15$. Of course this corresponds to a certain symmetry condition, but note that its periodic continuation is not even around $n=0$, as required by $(1)$ and $(2)$.

A simple example of a triangular signal similar to your example but with a purely real-valued DFT (i.e., satisfying $(2)$) is

x = [0 1 2 1]

Its purely real-valued DFT is given by

X = fft(x)
X =

   4  -2   0  -2

What follows is an explanation why the DFT coefficients $X[k]$ are real-valued if (and only if) condition $(1)$ is satisfied. If $X[k]$ is real-valued, we have

$$X[k]=X^*[k]\tag{4}$$

With

$$X[k]=\sum_nx[n]e^{-j2\pi nk/N}\tag{5}$$

where the sum is over one period, i.e., over $N$ adjacent indices, and where we assume that $x[n]$ is continued periodically, such that $x[n]=x[n+N]$. Using Eq. $(5)$ we can write $(4)$ as

$$\begin{align}\sum_nx[n]e^{-j2\pi nk/N}&=\sum_nx^*[n]e^{j2\pi nk/N}\\&=\sum_nx^*[N-n]e^{j2\pi (N-n)k/N}\\&=\sum_nx^*[N-n]e^{j2\pi k}e^{-j2\pi nk/N}\\&=\sum_nx^*[N-n]e^{-j2\pi nk/N}\tag{6}\end{align}$$

where I've used the index transformation $n\rightarrow N-n$, and the fact that $e^{j2\pi k}=1$. From $(6)$, Eq. $(1)$ follows immediately. Note that $(1)$ is dual to the fact that the DFT of a real-valued sequence is conjugate symmetric. In the case we've considered here, we have a conjugate symmetric sequence, leading to a real-valued DFT.


2D Convolution in MATLAB Causes Artifacts (Boundary Issues) I`m trying to do a 2D fast convolution in Matlab of large matrices. If I use FFT version based on convolution theorem ( https://en.wikipedia.org/wiki/Convolution_theorem ), there are some artefacts in the image. Only imfilter produces correct results, but it is working ~100x slower, than the FFT version (conv2fft).

What can be done to avoid these artefacts? I do know that in 1D we need to zero-pad, but here both the image and the kernel have the same size.

% generate image
len = 2^10;
CICcut = zeros (len);
CICcut = imnoise (CICcut, 'salt & pepper', 0.0001);
CICcut = CICcut.*(rand(len)).^2;
gauss = fspecial('gaussian', round(sqrt(len)), sqrt(sqrt(len)));
CICcut = imfilter (CICcut, gauss, 'replicate', 'conv');

% generate kernel
g = zeros(len);
lenMone = len-1;
for i = 1:len
    for j = 1:len
        g(i, j) = ((i-1)/lenMone - 0.5)^2 + ((j-1)/lenMone - 0.5)^2;
    end
end
g = -log(sqrt(g));

% convolution
tic
filtered    = imfilter (g, CICcut, 'replicate', 'conv');
toc
tic
filteredFFT = conv2fft (g, CICcut, 'same');
toc
tic
filteredN   = convn (g, CICcut, 'same');
toc

% display
figure('units', 'normalized', 'outerposition', [0 0.25 1 0.5])
subplot 151, imshow (CICcut, []); title ('Mass density')
subplot 152, imshow (g, []); title ('Green`s function')
subplot 153, imshow (filtered, []); title ({'Gravitational potential' 'imfilter'})
subplot 154, imshow (filteredFFT, []); title ({'Gravitational potential' 'conv2fft'})
subplot 155, imshow (filteredN, []); title ({'Gravitational potential' 'convn'})

Best regards,
Alex

If I understood right this is a duplicate question but I cannot find the link let me re-type the answer. Given that you want to implement a 2D discrete convolution of two images of sizes $N_1 \times N_2$ and $M_1 \times M_2$ by the method of DFT multiplication, then in order to avoid time (space) domain aliasing due to circular convolution implied by the inverse DFT stage, the sizes of the forward DFTs should be properly selected (in other words original time/space domain signals should be padded with zeros to increase thier sizes to overcome the aliasing).

Particularly the 2D forward and inverse DFT sizes should be selected as:
$$ L_1 \geq N_1 + M_1 -1$$ and $$ L_2 \geq N_2 + M_2 -1$$ in order to avoid circular artifacts and get the exact convolution.

Select the central portion of the resulting convolution to get the final image the same size of the original image.


Get position of different frequencies in sound signal Assumed I have a signal where different frequencies appear depending on the time position, i.e. a signal of $440\,\textrm{Hz}$ appears at $t=0$, a signal of $441\,\textrm{Hz}$ appears at $t=1$, a signal of $442\,\textrm{Hz}$ appears at $t=3$. Is it now possible to find out when the different signals appeared, and when they leave again, and relate that to the time position? The normal fft in one dimension just tells me that there are different frequencies, but not which.
My application is, that I have a sound file taken from an interferometer. The sound is generated by the interfering light beams. Different wavelengths generate different interference frequencies, but at different times. Thus I want to get the different frequencies in the resulting signal, and the time when they begin/end. The interference pattern generates a sine wave while interfering.

I don't understand the interferometer setup but have a look at:

Eric Jacobsen & Richard Lyons. DSP Tips & Tricks: The Sliding DFT. IEEE Signal Processing Magazine, March 2003, pp 74–80.

Eric Jacobsen & Richard Lyons. DSP Tips & Tricks: An Update to the Sliding DFT. IEEE Signal Processing Magazine, January 2004, pp 110–111.

You can use the sliding DFT (discrete Fourier transform) to calculate a windowed DFT starting at the time of each sample, not for all frequencies but for just the frequency you are interested about. See the subsections about generalized cosine windows in Wikipedia for the kind of windows that work with sliding DFT. Perhaps you can then use the magnitude and phase information from sliding DFT for your task.


What is the meaning of half window functions? In fig. 4 on p. 4 of the paper Extracting the frequencies of the pinna spectral notches in measured head related impulse responses (J. Acoust. Soc. Am. 118, 2005 ), a half Hann window is used (Hann window for positive time, zero for negative time). What is its meaning?

Because the head related impulse response (HRIR) is zero for negative time, you may as well interpret it so that they have used a full Hann window.


Proving DSP Sampling Theorem Schaum's Outline, Digital Signal Processing, Second edition, 2012, page 101:

Prove that:

$$X(e^{j\omega}) = \frac{1}{T_s}\sum_{k=-\infty}^{\infty}X(j\frac{\omega}{T_s}-j\frac{2\pi k}{T_s})$$

Solution

Impulse train, s(t), of period Ts is used to sample signal x(t):

$$s(t) = \sum_{n=\infty}^{\infty}\delta(t-nT_s)$$

And the sampled signal is:
$$x_s(t) = x(t)s(t)$$

Converting s(t) to frequency domain:
$$F\{s(t)\}= F\{\sum_{n=\infty}^{\infty}\delta(t-nT_s)\}$$

$$S(j\Omega)= \Omega_s \sum_{n=\infty}^{\infty}\delta(\Omega-n\Omega_s)$$

Where:

$$ \Omega_s = \frac{2\pi}{T_s} $$

$$ T_s = \frac{\Omega_s}{2\pi} $$

Now we convert the sampled signal to frequency domain:

$$
\begin{aligned}
F\{x_s(t)\} &= F\{\ x(t)\ s(t)\ \} \\
\\
X_s(j\Omega) &= \frac{1}{2\pi}\ \left[\ X(j\Omega)*S(j\Omega)\ \right] \\
\\
X_s(j\Omega) &= \frac{1}{2\pi}\ \left[X(j\Omega)*\Omega_s \sum_{n=-\infty}^{\infty}\delta(\Omega-n\Omega_s)\right] \\
\\
X_s(j\Omega) &=  \frac{\Omega_s}{2\pi}\ \left[X(j\Omega)* \sum_{n=-\infty}^{\infty}\delta(\Omega-n\Omega_s)\right] \\
\\
X_s(j\Omega) &= \frac{1}{T_s}\ X\left(j\Omega-jn \Omega_s\right) \\
\\
X_s(j\Omega) &= \frac{1}{T_s}\ X\left(j\Omega- \frac{j 2 \pi n}{T_s}   \right) \\
\end{aligned}
$$

Discrete-Time Fourier to Continuous-Time Fourier relationship:

$$x[n] = x(t=T_s n)$$

Take Fourier transform with respect to n to both side

$$ F_n \{ \text{ } x[n] \text{ } \} = F_n\{ \text{ } x(n T_s) \text{ } \} $$
$$X(e^{j\omega}) = \frac{1}{T_s} X\left(\frac{j\Omega}{T_s}\right) $$

Comparing $X_s(j\Omega)$ to $X(e^{j \omega })$, we see that $X_s(j \Omega)$ is a frequency-scaled version of $X(e^{j \omega})$ with the scaling defined by $\Omega=\omega/T_s$. The scaling that makes: $X(e^{j\omega})$ periodic with a period of: $2\pi$ is a consequence of the timing scaling that occurs when $x_s(t)$ is converted to $x[n]$

$$
X(e^{j \omega})=\left.X_s(j\Omega)\right|_{\Omega=\omega/T_s} = \frac{1}{T_s}\ X\left(\frac{j\omega}{T_s}-\frac{jn\omega_s}{T_s}\right)
$$

Properties Needed:

Property 1
$$F\{f_1(t)f_2(t)\} = \frac{1}{2\pi}[ F_1(j\Omega) * F_2(j\Omega)] $$

Property 2

$$F\left\{\sum_{n=-\infty}^{\infty} \delta(t-nT)\right\} = \Omega_s \sum_{n=-\infty}^{\infty}\delta(j\Omega - jn\Omega_s)$$
$$\text{   }Where: \Omega_s=\frac{2\pi}{T}$$

Property 3

$$a\text{ }(f_1(x)*f_2(x)\text{ } )=f_1(x)*af_2(x)=af_1(x)*f_2(x)$$

Property 4

$$f(x)*\delta(x-a) = f(x-a)$$

Property 5: scaling

$$ F\{ f(at) \} = \frac{1}{|a|} F(\frac{j\Omega}{a})$$

Graphs

# From Jupyter Notebook
%matplotlib inline
from numpy import *
from matplotlib.pyplot import *

def FourierSigCT(omega, omega_o=2.5):
    #   this Function "emulates" the spectrum of
    #   a fourier transformed CT signal with max frequency
    #   at +/- omega_o.
    #   Its created by taking a snipplet of a Cosine
    #   between -pi/2 and pi/2, and scaling it to
    #   the cutoff frequency in width..
    cutoff=(pi/2) / omega_o
    if ((omega < omega_o) and (omega > -omega_o)):
        return cos(omega*cutoff)
    else:
        return 0

def FourierSigSampled(omega_o, omega_s):
    t = arange(-10, 10, 0.01)
    n = arange(-5, 5, 1)
    Ts = 2*pi / omega_s
    y = zeros(len(t))
    for j in range(0,len(n)):
        for i in range(0,len(t)):
            y[i] = y[i] + (1/Ts)*FourierSigCT(t[i] - n[j]*omega_s, omega_o=omega_o)
    title("omega_o=%g  omega_s=%g" %(omega_o, omega_s))
    plot(t,y,'.')

t = arange(-3*pi, 3*pi, 0.1)
y = zeros(len(t))

for i in range(0,len(t)):
    y[i] = FourierSigCT(t[i])
title('X(j Omega)')
plot(t,y,'.')

FourierSigSampled(omega_o=2.5, omega_s=2.5*2)
FourierSigSampled(omega_o=2.5, omega_s=2.5*1.5)
FourierSigSampled(omega_o=2.5, omega_s=2.5*3)


Inverse DTFT of $e^{-j\omega}$ I've been trying to solve the following problem, but I noticed that I get two different answers by means of two different techniques:

Find the impulse response, $h(n)$, of the following system, described in frequency:
$Y(e^{j\omega}) = e^{-j\omega}X(e^{j\omega}) + \frac{dX(e^{j\omega})}{d\omega}$

Try #1: I found the difference equation for the system and iterated while making $x(n) = \delta(n)$:

$y(n) = x(n-1)+\frac{n}{j}x(n)$

$y(0) = 0 + 0$

$y(1) = 0 + 1$

$y(2) = 0 + 0$

And so on. Which yields $h(n) = \delta(n-1)$.

I found this same result by noticing that $X(e^{j\omega}) = 1$ for $x(n) = \delta(n)$, which gives: $Y(e^{j\omega}) = e^{-j\omega}$ and by using the time shift property we find that $h(n) = \delta(n-1)$.

Try #2: I used the inversion formula for $Y(e^{j\omega}) = e^{-j\omega}$:

$y(n) = \frac{1}{2\pi}\int_{-\pi}^{\pi}e^{-j\omega}e^{j\omega n}d\omega$ which yields $y(n) = \frac{sin(\pi(n-1))}{\pi(n-1)}$

Did I miss something? What is the correct result?

Your second try yielded a normalized sinc function. That function has a zero-crossing whenever $\pi(n-1)$ ia a nonzero multiple of $\pi$, which is for every $n\ne 1$ (assuming n is an integer) and is equal to 1 when $\pi(n-1)=0$.
Therefore, both solutions are the same.


Take FFT of audio represented in PDM I have a PDM data input from the MEMS microphone on the Pynq-Z1 development board. My eventual goal is to take the FFT of the audio samples and display the frequency bins on some LED.
I am struggling to understand what exactly is needed to convert the PDM signal to an acceptable input for the FFT. I have seen some videos explaining converting the PDM to PCM by applying a low pass FIR and then decimating the sampling rate.
Is this needed for my hobby project? Can I just average the PDM samples ( 0 to 1 ) over chunks of time T, to give a really simple idea of the amplitude at that time, and then feed these samples to the FFT?
How much do I need to process my PDM? I am very weak in the area of DSP.
Any help much appreciated.

In a PDM microphone, the sigma-delta convertor pushes the noise to frequency regions above the 0-20kHz spectrum. So if you'd take the FFT of the signal that comes directly from the PDM microphone, the straight zeros and ones, and only look at the bins that fall in the 0-20kHz region, you'd get the data that you need.
You can see this here, where a 16kHz sine wave has been converted to a PDM signal:

The noise floor on the left of the green 20kHz line is well below -70dB, but quickly increases to the right of it.
In the plot above, the PDM data rate is 768kHz. A common PDM microphone data rate will be more on the order of 3MHz. The problem with this is that you'll doing FFTs of a large amount of data and that you're throwing away a lot of unused data.
A moving average filter would indeed quickly reduce the number of data samples that go into the FFT. It will result in attenuation in the pass band that is too high for good audio quality, but it should be totally fine for an LED display.
I wrote about my own PDM microphone learnings in a series of blog posts.


What is the value of $0.01\log_{10}$? I have posted this question in math.stackexchange.com, [can be found here], which I originally thought it is a math question. I now believe this would fit more to the Signal Processing community. For convenience, I am reposting the question in below. I would appreciate any help.

  I am trying to understand what it means by $0.01\log_{10}$, which was published in the following paper at page 7 in IEEE.

  Berrier, Keith L., Danny C. Sorensen, and Dirar S. Khoury. "Solving the inverse problem of electrocardiography using a Duncan and Horn formulation of the Kalman filter." IEEE Transactions on biomedical engineering 51.3 (2004): 507-515.

  I cannot convert $0.01\log_{10}$ to a number since I expected to see something like $0.01\log_{10}^3$ to be able to say it is equal to 0.0048. The authors uses this expression twice in the paper.

  Could someone kindly help me in understanding what the authors meant by $0.01\log_{10}$?

  The following is a snapshot of the paragraph of the paper in which $0.01\log_{10}$ was used.

Differences in log scales correspond to ratios of the underlying values.  Saying you are within .01 on a $\log_{10}$ scale means the true value is within $10^{-.01}$ and $10^{.01}$ of the stated value, ratio wise.  That calculates out to 0.9772 and 1.0233.  Since the value is so close to 1, an interval on the ratio scale is close to an interval on a difference scale, so you could also say the results are within about 2.3%.

I haven't read the paper.

Hope this helps.

Ced

Followup:

So, on a whim I looked up "centibel" and found https://en.wiktionary.org/wiki/centibel.  A "cB" if you will.  Not too many other references though.

I should have also looked at the OP's math link.  The correct answer was already there.


How does minimum-latency partitioned convolution reverb work when you receive input samples in chunks, rather than one at a time? I'm writing a reverb system where I receive an input block of samples 480 elements long, do some operation on them, and pass the block on to the next effect.

I've been reading up on partitioned convolution reverb, and it seems to operate on the principle that you just accumulate samples until you have enough to convolve a short block of the impulse while the rest of the input comes in, getting bigger and bigger every chunk. How does this work if I can only start working on the input once it's all accumulated? How do convolution reverb plugins in DAWs work with this limitation (buffer sizes are usually lower, but never single-sample)?

For example, my impulse is 16384 samples long. I partition it into 480 + 32, then 480 + 64, 480 + 128, 480 + 256 (rounded up to 1024), etc. All the partitions add up to a length longer than it would have been if I'd just done 16384 + 480 samples (rounded up to 32768).

In my specific case, the source audio is a sample being played back at variable speed so I guess I could apply the reverb as each sample is being added to the buffer or even bake the reverb if I don't mind the room characteristics shifting with the pitch, but I'd like to understand how others do it first.

This is typically done using a segmented overlap add method or sometimes also refered to as a block convolver.

Let's assume your block size if 512 (makes the numbers a little easier).

Chop up your impulse response into 32 blocks of 512 samples each. Zero pad each block to 1024 samples and FFT. You know have 32 filters $H_0(z) ... H_{31}(z)$
On each new input block: zero pad to 1024 samples and FFT. Keep the most recent 32 signal spectra around, so you have $X_N(z), X_{N-1}, ... X_{N-31}(z)$
Multiply the spectra and sum them, $Y(z) = \sum_{k=0}^{31} H_k(z) \cdot X_{N-k}(z)$
Inverse FFT the result: you get 1024 time samples. The last 512 are the "overlap". Store for the next block.
Take the first 512 samples of your inverse FFT and add the overlap from the previous block. That's your output.

Repeat steps 2-5 for each new block.


How to measure acoustic pressure level under water with a hydrophone? this is my first post here and I am also new in signal processing.

Let's say I have a hydrophone with amplifier which is connected to a sound card plugged in to a computer. I gather the sounds to a WAV file every 3 seconds and analyze each file in a Python program. The system is supposed to detect loud events which are dangerous for sea animals.

I would like to know how to reliably measure acoustic pressure (loudness?) of events occuring under water based on WAV data. I am not sure if I should use exactly the same formulas as measuring SPL in air. I do not have the hydrophone yet so I can not really tell anything about it at this moment, but I am pretty sure the catalog notes of the hydrophone will be available for me when I finally get it.

It would be good enough if the measurement had an error under 5 dB SPL. It does not have to be perfect, I just want it to be relative and adapted to under water characteristics.

Could the simple measurement be calibrated by shifting measures by single constant in dB?

I found a similar solution in Determine Sound Intensity in SPL. But the problem is I do not know if I will be able to have all the data needed for this, mainly about amplification levels. Would the calibration-like solution be good in this case? How can I do it?

A hydrophone should come with at least a typical response curve.

The thing to aware of is that that response represents an idealized far field response in an infinite space referenced to an wave of an ideal point source taken at one meter (or yard).

A good reference that covers the conversion of acoustic to electrical energy is

  Horton, Joseph Warren. Fundamentals of sonar. United States Naval
  Institute, 1959.

but his book is scarce as hen's teeth and requires luck or access to a good University library that collected physics and engineering books from the 1950's.

There are other books on transducers.  Essentially an accurate SPL measurement is a bit more than a hydrophone hooked into a sound card.  The field expedient calibration  method is to use a calibrated source at a known direct path distance at least 100 wavelengths in distance.

In most cases the cable you need is going to be too long for a unamplified transduction signal.  It's best to have the converter at the sensor and send over a digital signal but have weak analog signals close to digital circuitry.

Long cables usually imply higher voltages to minimize resistance losses.  Sharks are attracted to cables that have high voltages so don't be surprised when you see cables offering Kevlar with steel reinforcement.

Ideally access to pressure chambers to test your the water-tight integrity of your hyrophone cable is a very good thing to have if you expect to dangle more than 20 feet and the deeper you go the more important it is to pressure check your set up.


What Does It Mean Exactly When Two Parts of a Signal Are Correlated? I stumble quite often upon the notion that two or more parts of a signal are correlated to describe semi-formally that they belong together. For example in image processing, two pixels on an edge feature tend to be correlated whereas two adjacent parts of a 3D structure that represents water droplets in a particle simulation are less correlated. My question is what is the exact idea behind this notion.

Yeah, it can mess you up pretty badly if you don't get the fundamentals right off the get-go.  This is how I interpret correlation, and it has worked for me for what I do for a living.

Let's start off with a relatively simple example.  Take a look at the following figure (pulled from dspguide... this is actually a great online book for knowing the basics of DSP).

We have an antenna that transmits a short burst of radio wave energy in some direction. If the propagating wave strikes an object.... like a helicopter in this figure, a small fraction of the energy is reflected back toward a radio receiver.  This receiver is close to the transmitting antenna.

This short burst of radio energy, for the sake of this example, is a small triangular shape.  When the signal is reflected off of the helicopter, and then echoed back to the receiver, this signal will consist of two parts:

A shifted and scaled version of the transmitted pulse, and
Random noise, resulting from interfering radio waves, thermal noise in the electronics and other factors.

Loosely speaking, we can actually figure out how far the object is by using this concept.   Since radio signals travel roughly at the speed of light, the shift between the transmitted and received pulse is a rough measure of the distance to the object being detected.

As such, this is our general problem:

  Given a signal of some known shape, what is the best way to determine where (or if) the signal occurs in another signal?

The best way to answer this is correlation.

There are two different paradigms for computing correlation.  The first one is called auto-correlation, where you are comparing a signal with shifted time offsets of itself.  This paradigm that we are describing (also seen in the figure) is defined as cross-correlation, where we are comparing with another signal, notably the received signal.  We essentially are comparing the received signal with shifted versions of the original transmitted signal.  Basically, we take a look at what we have received and what was transmitted.  We take what was received, and time shift the original transmitted signal over by different time values.  We then do a comparison with each of these signals and the received result.  Whichever gives us the highest value will denote how far away the helicopter is.

The amplitude of each sample in the cross-correlation signal is a measure of how much the received signal resembles the target signal, at that location. This means that a peak will occur in the cross-correlation signal for every target signal that is present in the received signal. In other words, the value of the cross-correlation is maximized when the target signal is aligned with the same features in the received signal.

If there is noise on the received signal, there will also be noise on the cross-correlation signal. It is an unavoidable fact that random noise looks a certain amount like any target signal you can choose. The noise on the cross-correlation signal is simply measuring this similarity. Except for this noise, the peak generated in the cross-correlation signal is symmetrical between its left and right. This is true even if the target signal isn't symmetrical.

A good thing to remember is that the cross-correlation is trying to detect the target signal, not recreate it. There is no reason to expect that the peak will even look like the target signal.  Correlation is the optimal technique for detecting a known waveform in random noise. To be perfectly correct, it is only optimal for random white noise. Using correlation to detect a known waveform is frequently called matched filtering.

tl;dr - Correlation is a measure of how much one signal resembles another.  The signal can be images, features, edges, etc.  It is simply a measure of resemblance between one signal and another.


Draw outline of a text in an image I need to draw an outline of text in an image. Consider the below image as an example.

Can anyone suggest me a direction to achieve this? I don't have experience in image processing. Any tools/lib or anything that can generate this kind of image will be helpful.

You just need to binarize the image then find the hull. In Mathematica it looks like this:

<< ComputationalGeometry`
img = EdgeDetect@Rasterize["managers. In\n more lean\n result of\nss\nkforce and\nt \
points and"]
pts = Position[ImageData@img, 1]
pgp = PlanarGraphPlot[pts, ConvexHull@pts][[1, 2]]
Graphics[{Black, Point@pts, pgp}]

If this is inadequate I would find the convex hull of the text two lines at a time, then fuse the results.


Why is $\int^\infty _{0^-}\delta(t-nT)e^{-st}dt = e^{-nsT}$? I'm currently in the process of going over the $\mathcal Z$-transform and more specifically its derivation. I understand and I am able to follow it up until the final step whereby involving the Laplace transform.

The step that is confusing me is indicated as the step within the red box, derivation taken from Wikipedia page on the $\mathcal Z$-transform. I do not understand how this integral solves to $e^{-nsT}$. I have searched for possible identities relating to this but have not been able to find any help.

\begin{align}
X_q(s)&=\int_{0^-}^\infty x_q(t)e^{-st} dt\\
&=\int_{0^-}^\infty \sum_{n=0}^\infty x[n]\delta (t-nT)e^{-st} dt\\
&=\sum_{n=0}^\infty x[n]\color{red}{\boxed{\color{black}{\int_{0^-}^\infty \delta (t-nT)e^{-st} dt}}}\\
&= \sum_{n=0}^\infty x[n]\color{red}{\boxed{\color{black}{e^{-nsT}}}}
\end{align}

Thanks for taking the time for looking at my question. Any help would be greatly appreciated!

Regards

Remember the definition of the delta function: it integrates $1$ all over the t-axis but it is zero for any $t\neq 0$. This means that if we multiply any function by $\delta (t)$, then the function is going to be zero everywhere except for $t=0$. In this case, we have

$$ \delta (t-nT)e^{-st}=
\left\{
 \begin{array}{ll}
  \delta (t-nT)e^{-snT}  & \mbox{if } t = nT \\
  0 & \mbox{if } t\neq nT
 \end{array}
\right.$$

So if we integrate this function, we have that the integral is zero for every point except for $t=nT$. So

$$\int_0^\infty \delta (t-nT)e^{-st}dt=\int_0^\infty \delta (t-nT)e^{-snT}dt=e^{-snT}\int_0^\infty \delta (t-nT)dt$$

But by definition, $$\int_0^\infty \delta (t-nT)dt=\int_{-\infty}^\infty \delta (t-nT)dt=1$$

So $$\int_0^\infty \delta (t-nT)e^{-st}dt=e^{-snT}$$

Sidenote: For this to be true, $n$ must be positive. If it is negative, then the integral will be $0$ (as the function inside the integral will be zero for all the integration interval).


Does Scrambling increase speed of signal transmission overall? And is Scrambling a redundant process? Title is a badly worded question, on which I do not have enough knowledge on the subject to answer on my own.
Anyone to help out?

Scrambler is a randomizer with a purpose to encode transmitted data so that there will be no long sequences of only zeros or only ones in transmitted data. It is done to support proper timing recovery, when receiver waits for bit transitions to synchronize. Scrambling is usually done for channels with self-synchronized line codes, e.g. ISDN u-interface and Ethernet.
Generally, scrambler doesn't need to change the amount of transmitted bits, hence it doesn't change the speed of data transmission.
You can compare a scrambler to error-detection and correction schemes, which add redundant bits to support its functions. You can similarly compare it to 8b/10b encoding, when 2 bits are added to every 8 bits of data. As a scrambler doesn't add additional information to the data, it is not a redundant process in that sense.
Also, probably irrelevant in most cases, but it can be considered redundant, when used for a signal that inherently doesn't need scrambling. For example, UART uses short transactions and receiver is capable to get all bits only based on knowledge of a baud rate, so use of a scrambler would be redundant.


How to compute two Dimensional SNR The signal is the white dot or it is a point target. and I have to find the SNR of this image/point target. As i increase noise in thie image, black back ground shows some disturbance. AnHow should i approach this problem. I have to do it in matlab code.

Basically I am working on the compressed response(poin target) of an image. Here in image i have a point target(white bright dot) and in zoom your can see its RGB cells. I want to compute 2 Dimensional SNR of point target. Point target is the signal and rest of it is noise.
This processed image has an array of (512x1024)rows*coulums= 524288 pixels
My approach is follow the basic formula SNR= 10log(Psignal/Pnoise)
The point target (white bright point) has a main lobe(compressed response) (bump)
The rest are side lobes.
##Power_signal= Main_loob Energy/pixels covered by point target##Power_noise= Rest side lobes/(Total Pixels - pixels of point target)##
I am working on this approach and not completed yet...
I need guidance on 2D SNR:

Is my approach correct
2)Do we consider back ground as noise or not
3)Any other methods that will make the problem easy to solve


Find filter coefficients to model a device using its measurement I am trying to realize a digital filter that has the same freq. response of an existing speaker. I have fed an audio sine sweep to the speaker and measured the speaker output, both at 48kHz. Then I perform FFTs of the input(X) and output(Y), divide the absolute values point-wise (absY./absX) to get the transfer function/freq. response (H).

Now I would like to determine the filter coefficients B and A for an FIR/IIR filter, so that I can model the speaker's response digitally. I understand my transfer function is in the domain of w, and filter coefficients are the coefficients of the difference equation of the filter.

This is where all my text-book theory seems to fall apart, and I was hoping someone could clarify:

If I understand correctly, since my data is discrete (sampled at 48k), my data is in terms of 'n'. If I performed an ifft of H (which is complex), the resulting vector is basically my filter coefficients?
Matlab suggests using invfreqz to find filter coefficients, given a complex frequency response. Does this function map from w to z domain? If so, am I right in understanding that the FFT needs to be converted to Z domain to derive its filter coefficients?
What parameters does a auto-regressive filter fit to? for eg: Yule_Walker, or Levinson, or if I want to run a gradient descent to fit a filter? What is the error calculated between?

We can represent a digital filter as:

A difference equation, relating input and output at every time step $n$:
$y[n] = \sum_{k=0}^{K-1} b[k] x[n-k] - \sum_{m=1}^{M} a[m] y[n-m]$. The filter coefficients are $b_k$ and $a_m$. In case of an FIR, $a[m] = 0\,\forall\, m $. This form is basically how the filter can be implemented (in the discrete time domain).
By z-Transform of the difference equation, we can also represent the above equation as: $Y(z) = \frac{\sum_{k=0}^{K-1} b[k] z^{-k}}{1 + \sum_{m=1}^{M} a[m] z^{-m}} X(z) = H(z)X(z)$.
The frequency response is simply this z-Transform evaluated at the unit circle $H(z = \mathrm{e}^{j\omega})$ (this is often simply written as  $H(\omega)$).
The impulse response is the inverse transform of $H(\mathrm{e}^{j\omega})$. By convolution with the impulse response, we also can get the output of the filter, as in $y[n] = h[n] * x[n]$. However, $h[n]$ is only finite for an FIR (hence the name, obviously), where we can easily check that in fact it is simply the coefficients $b[k]$.

I hope this helps you with understanding the big picture. So to answer your questions:

The IFFT of the frequency response yields the impulse response. This corresponds directly to the filter coefficients only in the case of an FIR filter. In case you want to model your system as an IIR, you can use an algorithm such as the kind implemented in MATLAB's invfreqz, which is essentially an optimization to fit your frequency response to the complex rational function $H(z)$.
Yes, since $H(z)$ is a representation of the filter's difference equation, the coefficients of which are what we call filter coefficients.

In any case, these are important basics you can learn from e.g. the classic book on DSP. Even the MATLAB documentation (which im sure you have already been checking) is a decent resource.


Does Wav format supports writing complex values I have a sine tone, doing some math, and finally gets an harmonically accumulated sine tone in complex numbers. I do want to preserve the phase information of imaginary part on the other end, which will continue processing the complex signal. So,

Is there any way i can safely pack the complex number into any of standard file formats, such as wav ?

Note : Currently I am  writing this signal's byte values to a .bin file

You could create a standard 2 channel wav file. 1 channel for the real part and 1 for the imag. Wav files supports both 16bit integer, 32bit integer and 32bit floats etc. The channel data will be interleaved (I believed it's called). Meaning that your data will be written: [REAL][IMAG][REAL][IMAGE]...
Prolly just like your own bin file.

Standard programs like Audacity will be able to read the file and make calculations on each (or both) channels.


FFT Plot Axis Confirmation I have a 4096 sequence of FFT coefficient. All of which a real values. My sampling frequency is also 4096 samples per seconds. Now I want to make a plot. I know my Y axis would be my FFT magnitude values, but what would my x- axis be. Would it be my sampling frequency, that is 1,2,3,4.........4096?

Here's a short snippet of code that returns the domain of the FFT depending on if you have the onesided or twosided FFT. This is from some code I posted on a somewhat related question here.

% Nfft - FFT size
% Fs - Sampling frequency in Hz
% oneside - flag indicating FFT results will be onesided or two sided

% Compute the frequency vector
if oneside == true
  f = fs * (0:Nfft/2).'./Nfft;
else
  f = fs * (-Nfft/2:Nfft/2-1).'./Nfft;
end


low frequency transform I was wondering which type of signals have bounded-support Fourier transform. e.g. their transform is limited from zero to some non-infinite frequency.

The main reason I'm into it is that in Shannon's sampling theorem, we know that for some sampling rate, reconstruction is possible without aliasing, but I guess that's only true for functions which have a limited (cut-off) frequency

Only waveforms that have infinite support in one FT domain have finite support in the other.

This implies that only impossible (in the real world) or fictitious signals comply perfectly with the sampling theorem (at a finite sampling rate).  However given a non-zero noise floor (finite sized samples, quantum/Planck sampling time jitter, etc.), a common assumption is that the portion of the infinite spectrum sufficiently around or below the noise floor can be treated as zero, thus one can pretend that these signals are band-limited for practical purposes.


Transfer function of linear phase filter & its frequency domain representation Text of my exercise request:

  Determine the transfer function $H(z)$ of a causal linear phase FIR filter with zeros at $z= \frac{1}{3}$ and $z=-2$. The value of the impulse response at $n=0$ equals $2$.

So I calculated $H ( z )$, taking into account the zeros and their inverses (due to the linear phase property): $$H(z)= k(z^{-1}-3) (z^{-1}-\frac {1}{3})(z^{-1}+ \frac{1}2)(z^{-1}+2)$$
intermediate steps of the multiplication, take me
$$  H(z) = k (z^{-4} - \frac{5}{6}z^{-3} - \frac{19}{3}z^{-2} -\frac{5}{6}z^{-1}+1) $$
and finally $$ h(0) = k=2 $$

but what is the corresponding DFT? I think it's a 5-point DFT, but I do not know calculate it, I just say that
$$H(K)=H(z) _{|z= e^{\frac{j2 \pi k }{5}}}$$
but it would be necessary to do the calculations?

You finally have changed your $H(z)$ and added a qualifier "linear phase" , so is the new answer:

Now assuming that you have already determined $H(z)$ from the problem specifications, so the rest of obtaining DFT $H[k]$ can be approached in a number of ways one of which (perhaps longer than necessary) is the following:

Using the linearity and shift properties of the Z-Transform, we can show that that if a given signal $x[n]$ is of the form $$x[n] = a\delta[n] + b\delta[n-1] + c\delta[n-k]$$ then its Z-Transform is of the form: $$X(z) = a + b z^{-1} + c z^{-k}$$

Applying this observation, in the reverse direction, to $H(z)$ in your problem we get a prototype for $h[n]$ as: $$h[n]=k(\delta[n] - {5\over6} \delta[n-1] - {19\over3} \delta[n-2] - {5\over6} \delta[n-3] + \delta[n-4])$$

Also as you state that $h[0]=2$, yields $k=2$, hence the signal $h[n]$ is: $$h[n]= 2\delta[n] - {5\over3} \delta[n-1] - {38\over3} \delta[n-2] - {5\over3} \delta[n-1] + 2\delta[n-4]$$

Obtaining DFT $H[k]$ of this $h[n]$ is actually a very similar procedure to obtaining $H(z)$, you shall simply plug $e^{j{2\pi \over N} k}$ in place of $z$ where you simply get this:

$$H[k]= 2 - {5\over3} e^{-j{2\pi \over N} k} - {38\over3} e^{-j{4\pi \over N} k} -{5\over3} e^{-j{6\pi \over N} k} + 2e^{-j{8\pi \over N} k}$$

where $k$ ranges from $0$ to $N-1$ and plugging $N=5$ produces:
$$H[k]= 2 - {5\over3} e^{-j{2\pi \over 5} k} - {38\over3} e^{-j{4\pi \over 5} k} -{5\over3} e^{-j{6\pi \over 5} k} + 2e^{-j{8\pi \over 5} k}$$

Note that eventhough it is perfectly valid to state that, under suitable conditions, $H[k]=H(z)$ evaluated at $z=e^{j{2\pi\over N} k}$ , I guess it does not make an answer...


Relationship between sampling and number of pixels? In case of digital images, what is the relationship between sampling and number of pixels? For upsampling, the number of pixels increass, and for down sampling, the number of pixels decrease??

Yes, that's correct. It's a natural extension of one dimensional up and down Sampling. However, one can define lattices (arbitary shapes) for sampling images, so one could downsample more in a particular axis than the other. Essentially instead of a sampling period, we have a sampling grid. There are benefits of using different freqeuency grids. Of course all of this holds without aliasing only when the sampling theorem holds. Otherwise you get aliasing.


How do you deconvolve irregularly sampled points? $f:\mathbb{R}^n\rightarrow\mathbb{R}$ is a signal and $g:\mathbb{R}^n\rightarrow\mathbb{R}$ is a known point-spread function (say, a Gaussian).

A system samples $f\star g$ at a known sequence of irregular points, $(\mathbf{x}_n)_n$. Are there any methods, analytic or numerical, to deconvolve the samples from $g$ to reveal samples of $f$?

If $f$ were band-limited and the sequence of sample points $(\mathbf{x}_n)_n$ were periodic and finer than the sampling theorem's frequency, then I believe by the convolution theorem sampled-$g$'s  Fourier transform could be divided out pointwise from that of sampled-$f\ast g$, yielding sampled-$f$'s Fourier transform.

You can use a graph-based Signal Processing framework to apply regular deconvolution algorithms. The whole point is that traditional operators (such as FFT) can be defined on graph structures, and graph structures do handle "by design" without any overhead irregular sampling patterns.

Some useful reference scan be found on my former office mate's webpage.


Fixed-point Exponential Moving Average filter with given cutoff frequency I want to implement in fixed point arithmetic an exponential moving average filter, with a specific cutoff frequency for a given sampling rate.
The formula for the filter is
$$
y_n = \alpha x_n + (1 - \alpha)  y_{n-1}
$$
Looking at Exponential moving average cut-off frequency one can derive the value for alpha for a given sampling rate and a wanted cutoff. For instance, for a sampling rate of 1000 Hz and a cutoff of 50 Hz, the value of alpha is 0.267730531659313:

Fs = 1000
f3db = 50 %This is the cutoff frequency that I want
format long;
omega3db = f3db * pi/(Fs/2)
alpha = cos(omega3db) - 1 + sqrt(cos(omega3db).^2 - 4*cos(omega3db) + 3)  =  0.267730531659313

Now I tried to implement the exponential moving average filter in fixed point arithmetic as described here:  https://www.embeddedrelated.com/showarticle/779.php but I am getting an offset that I can't explain.

Here my Arduino code:

#define PERIOD_MICROSECS 1000

static unsigned long lastRead = 0;

int analog_pin = 0;   // sensor connected to analog pin 0
int analog_input0 = 0;  // variable to store the read value;
int analog_input0_filtered_fixed_point = 0;

/*
 * Formula to compute the alpha parameter of the EMA given a desired cutoff frequency.
 * Taken from: https://dsp.stackexchange.com/questions/40462/exponential-moving-average-cut-off-frequency
 *
 * Fs = 1000
 * f3db = 50 %This is the cutoff frequency that I want
 * format long;
 * omega3db = f3db * pi/(Fs/2)
 * alpha = cos(omega3db) - 1 + sqrt(cos(omega3db).^2 - 4*cos(omega3db) + 3)  =  0.267730531659313
 *
*/

//I calculate in advance the coefficients for the filter
int Q = 12; // I can choose Q freely
//I want a cutoff of 50Hz at a sampling frequency of 1000 Hz. So I calculate the corresponding alpha with formula above
//For cutoff 50 Hz, alpha = 0.267730531659313; then alpha_scaled = round(2^Q*alpha) = round(2^12* 0.267730531659313) = round(4096* 0.267730531659313) = 1097
int alpha_scaled = 1097;
// N ~= -log_2(alpha) = -log_2(0.267730531659313) = 1.901146423365515 ~ 2
int N = 2;

int low_pass_EMA_fixed_point(int x, int y, int alpha_scaled, int Q, int N){
  y += (alpha_scaled * (x-y)) >> (Q-N);
  return y >> N;
}

// ************************* //

void setup() {
  Serial.begin(115200);
  analog_input0_filtered_fixed_point = analogRead(analog_pin); //set EMA y for t=1
}

void loop() {
  if (micros() - lastRead >= PERIOD_MICROSECS) {
        lastRead += PERIOD_MICROSECS;

        analog_input0 = analogRead(analog_pin);
        analog_input0_filtered_fixed_point = low_pass_EMA_fixed_point(analog_input0, analog_input0_filtered_fixed_point, alpha_scaled, Q, N);

        //Check the original and filtered signals with the serial plotter
        Serial.print(analog_input0);
        Serial.print(" ");
        Serial.println(analog_input0_filtered_fixed_point);
  }
}

See above a plot with the serial plotter.

Any idea about how to solve the offset problem? What am I doing wrong? Any suggestion about how to solve the general problem of a fixed point exponential moving average filter with a specific cutoff frequency?

You should scale your analog input to match your desired resolution, it's better to have bigger resolution (specially when dealing with IIR filters), I don't know the reason of returning y>>N, probably the scaling problem is there, I would implemented it like this:

int low_pass_EMA_fixed_point(int x, int y, int alpha_scaled, int Q, int N)
{
  int32_t y_temp;
  y_temp = (alpha_scaled*(x<<N))+((1<<Q)-alpha_scaled)*y;
  return (int)(y_temp >> Q);
}

Also a good practice is test your filter with an impulse to check if it was implemented correctly before using the actual data.

EDIT: tested this program with onlinegdb.com (the printf are used to copy the result into octave):

#include <stdio.h>
#include <stdlib.h>

#define NUM 1024

int low_pass_EMA_fixed_point(int x, int y, int alpha_scaled, int Q, int N);

int main()
{
    int alpha_scaled = 1097; //Q12
    int Q = 12;
    int N = 2;
    static int x[NUM];
    int y = 0;

    x[0] = 1024; //input in Q10
    printf("y=[");
    for(int i=0;i<NUM;i++)
    {
        y = low_pass_EMA_fixed_point(x[i], y, alpha_scaled, Q, N);
        if(i == NUM-1)
        {
            printf("%d];",y); //print the last output
        }
        else
        {
            printf("%d,",y);
        }
    }
    return 0;
}

int low_pass_EMA_fixed_point(int x, int y, int alpha_scaled, int Q, int N)
{
  int32_t y_temp;
  y_temp = (alpha_scaled*(x<<N))+(((1<<Q)-alpha_scaled)*y);

  return (int)(y_temp>>Q);
}

And I'm getting this frequency response in octave:


Fourier transform of $\sum_{n=-\infty}^\infty(-1)^n\delta(t-nT_0)$ Given $x(t)$ and $h(t)=\sum_{n=-\infty}^\infty(-1)^n\delta(t-nT_0)$, I have to compute $Y(f)$, where $y(t)=x(t)h(t)$. I have thought about using that, in this case, $Y(f)=X(f)*H(f)$. I know that $\mathscr{F}(\sum_{n=-\infty}^\infty\delta(t-nT_0))=T_0^{-1}\sum_{n=-\infty}^\infty\delta(t-nf_0)$, but how can I deal with that $(-1)^n?$

HINT:

Note that the given $h(t)$ can be written as

$$h(t)=g(t)-g(t-T_0)\tag{1}$$

with some $g(t)$ the Fourier transform $G(f)$ of which you know. So from $(1)$ you then get

$$H(f)=G(f)\left(1-e^{-j2\pi fT_0}\right)\tag{2}$$


How to find Leq for a WAV audio file like an SPL Meter would? Given an audio file (WAV), I need to find the $\rm L_{A_{eq}}$ (continuous A-weighted equivalent sound pressure level) of the audio file. For this purpose, I am following this workflow:

Get a frame of audio (say $0.125$ seconds of samples)
Scale the samples to find the $\rm SPL \ re \ 20 \ \mu Pa$ based on calibration factor obtained from a calibration tone.
Apply A weighting to the frame of audio data.
For each frame of audio data, find the SPL (dB) as:

$$
\begin{align}
\textrm{pressure}_{\rm Ref}    &= 20\times10^{-6};\\
\textrm{amplitude}_{\rm rms}  &= \sqrt{\rm mean( audData_1frame.^2 )};\\
\textrm{dBspl}_{\rm perFrame} &= 20\times \log_{10}\left(\frac{\rm amplitude_{rms}}{\rm pressure_{Ref }}\right);\\
\end{align}
$$

This gives the $\rm SPL \ re \ 20 \ \mu Pa$ of a frame of data. This would be the Short $\rm L_{eq}$.

How do I proceed further to get these 3 quantities:

$\rm L_{eq}$ = The total continuous equivalent sound pressure level. This wold be a scalar value for the entire audio file. If it has to be updated in a plot, how do we integrate the SPL of each frame with the previous frames?
$\rm L_{eq_{fast}}$ = Sound pressure level with $125\ ms$ time weighting
$\rm L_{eq_{slow}}$ = SPL with $1\ s$ time weighting

how do we integrate the SPL of each frame with the previous frames?

You're measuring the RMS value of the (filtered) samples, which is sqrt(average(samples^2)), so if you're finding the RMS value in multiple chunks, it would just be sqrt(weighted_average(chunk1^2, chunk2^2, chunk3^2, ...)) where the average is weighted by the number of samples in each chunk (typically they are all the same length, but maybe a recording is not an exact multiple of chunk size, for instance).

To find the time weighting, the signal flow is:

Apply A-weighting filter
Square the samples
Apply a low-pass filter with a real pole at 1/τ (where τ is 125 ms or 1 s in your case)
Square root
20*log10() of the signal to convert to dB

It's the same as the RMS of the weighted signal, above, except that instead of an unweighted average over the whole signal, you're using a sliding-time average (same as https://en.wikipedia.org/wiki/Moving_average#Exponential_moving_average I think).

The time-weighted SPL will be a signal, not a scalar for each frame.


Concatenated Deletion & Sticky Channel Capacity Suppose you are sending information through a deletion channel with a deletion rate of .5 and then the output of this channel is sent through a sticky channel with a repeat rate of .5 also.  No error correction is allowed in this scenario, the information is not coded in any way.I am not looking for an error correction code for this composite channel, I am NOT sending any information through the channel, I just want to calculate the channel capacity. It is well known that the capacity of the IID deletion channel is unknown. What I am hoping for is to show that if you concatenate a deletion channel with a sticky channel that the capacity for the combined channels is zero, or close to zero. Anybody out there who can do this?

The identical and identically distributed (IID) binary deletion and sticky channels are both general repeat channels each characterized by a discrete probability distribution of repeating an input bit in the output $D$ times with $D \in 0, 1, 2, 3, 4, \dots$ with $D\!=\!0$ meaning that the input bit does not appear in the output of the channel at all. The question's probabilities for the deletion channel (first channel, hence the subscript $1$) are given by:

$$P(D_1) = \left[\tfrac{1}{2}, \tfrac{1}{2}, 0, 0, 0, \dots\right],\tag{1}$$

which is to say $P(D_1\!=\!0)=P(D_1\!=\!1)=\frac{1}{2},$ as given in the question.

An IID sticky channel with input duplication probability $\frac{1}{2},$ as given in the question, and of a variety that may at most double the number of bits, is also a general repeat channel, with $D_2\text{-}$times repeat probabilities:

$$P(D_2) = \left[0, \tfrac{1}{2}, \tfrac{1}{2}, 0, 0, \dots\right]\tag{2}$$

Generally, without hard-coding the above numbers, the concatenation (cascade, in series) of the above two channels into a single channel is also a general repeat channel with $D_1D_2\text{-}$times repeat probabilities:

$$P(D_1D_2) = \left[P(D_1\!=\!0) + P(D_2\!=\!0)\\- P(D_1\!=\!0)P(D_2\!=\!0),\\
P(D_1\!=\!1)P(D_2\!=\!1),\\
P(D_1\!=\!2)P(D_2\!=\!1) + P(D_1\!=\!1)P(D_2\!=\!2),\\
P(D_1\!=\!3)P(D_2\!=\!1) + P(D_1\!=\!1)P(D_2\!=\!3),\\
P(D_1\!=\!4)P(D_2\!=\!1)\!+\!P(D_1\!=\!1)P(D_2\!=\!4)\!\\+\!P(D_1\!=\!2)P(D_2\!=\!2),\\
\dots\right]\tag{3}$$

These are probabilities for the possible values of the product of independent variables $D_1$ and $D_2.$ For example the last entry before the ellipsis collects probabilities for the three mutually exclusive events $4\times1=$ $1\times4=$ $2\times2$ $=4.$
Note that $P(D_1D_2) = P(D_2D_1),$ so it does not matter which channel comes first. Using the values of $P(D_1)$ and $P(D_2)$ from Eqs. 1 & 2, we get:

$$P(D_1D_2) = \left[\tfrac{1}{2}, \tfrac{1}{4}, \tfrac{1}{4}, 0, 0, \dots\right]\tag{4}$$

This characterization of the concatenated channel appears more complex than the same for the deletion and sticky channels. For each of the latter, finding the capacity appears to be an open research problem (now at the end of 2018). So I doubt the capacity of the concatenated channel is known either. However, it may be possible to find upper and/or lower bounds. For the capacity of the deletion channel, there seems to be a lower bound $1/18,$ see: Capacity of a Binary Deletion Channel.

Your question says that there should be no error correction. As far as I understand, the channel capacity does not depend on error correction, but if a working code can be demonstrated, then it gives a lower bound for channel capacity.

Mahdi Cheraghchi calls the general repeat channels $\mathcal{D}\text{-}$channels, where $\mathcal{D}$ would be the distribution of $D=D_1D_2$ for the concatenated channel. Perhaps you will find something of interest in their work (Cheraghchi 2018):

  A $\mathcal{D}\text{-}$repeat channel, given a bit,
  draws an independent sample $D ≥ 0$ from $\mathcal{D}$
  and outputs $D$ copies of the received bit.

  Our framework can be applied to obtain capacity upper
  bounds for any repetition distribution (the deletion and
  Poisson-repeat channels corresponding to the special cases
  of Bernoulli and Poisson distributions).

References

Mahdi Cheraghchi. 2018. Capacity upper bounds for deletion-type channels. In Proceedings of the 50th Annual ACM SIGACT Symposium on Theory of Computing (STOC 2018). ACM, New D_2ork, ND_2, USA, 493-506. DOI: https://doi.org/10.1145/3188745.3188768 (arD_1iv:1711.01630 [cs.IT])


Deriving expression for the DTFT of a rectangular window Looking at the picture above, how did the author get from point A) to B)?

My Approach: Multiply A) by $e^{j\omega/2}/e^{j\omega/2}$. Now I am stuck with simplying the numerator.

There's a mistake in the denominator of B. This is the way it should be done:

$$\begin{align}R(\omega)&=\frac{1-e^{-j\omega (M+1)}}{1-e^{-j\omega}}\\&=\frac{e^{-j\omega(M+1)/2}}{e^{-j\omega/2}}\cdot\frac{e^{j\omega(M+1)/2}-e^{-j\omega(M+1)/2}}{e^{j\omega/2}-e^{-j\omega/2}}\\&=e^{-j\omega M/2}\frac{\sin\left(\frac{M+1}{2}\omega\right)}{\sin(\omega/2)}\end{align}$$


How to efficiently control an FIR's magnitude response by altering its phase spectrum Question:

Extensively searching the space of all possible vectors of length $n$ to satisfy a (non-overdetermined) requirement is possible in principle. Hence, there is a way to calculate a complex FIR kernel of length $n$ whose magnitude response goes through $2n$ given points, when its phase response is allowed to be completely arbitrary. (Corollary, a real FIR kernel with $n$ given magnitude points and $n$ mirrored points can be calculated.) The extensive search method is hopelessly inefficient, though. Is there a practical way to calculate an FIR with such requirements?

Context:

It is common DSP practice to use the DFT to convert an array of real magnitude points $M = \{m_0, m_1, \dots, m_{n-2}, m_{n-1}\}$, $m_i \in \mathbb{R}$ along with an array of real phase points $P = \{p_0, p_1, \dots, p_{n-2}, p_{n-1}\}$, $p_i \in \mathbb{R}$ to an FIR kernel $K = \{k_0, k_1, \dots, k_{n-2}, k_{n-1}\}$. $K$'s elements are generally complex, but if $M$ is symmetric (i.e. $m_i = m_{n-i}$ for $0 < i < \frac{n}{2}$) and $P$ is conjugate symmetric (i.e. $p_i = - p_{n-i}$ for $0 < i < \frac{n}{2}$), then $P$'s elements are completely real.

The magnitude response (absolute value of the DFT) of $K$ will always match $M$ at the individual points $m_i$. Between these fixed points, however, it can attain very unintuitive and undesirable values. This undesirable behavior can easily be visualized by appending a large number of zeros to $K$ and then calculating the absolute value of the DFT. Changing any $p_i$ changes the behavior of the magnitude response of $K$ in between the $m_i$. Thus, two FIRs derived from the same magnitude spectrum $M$, but with different phase spectra $P$, can have very dissimilar magnitude responses. This is illustrated by the figure below:

The upper plot shows three randomly chosen examples of $P$, called $P_0$, $P_1$ and $P_2$. The lower plot shows a randomly chosen example of $M$ and the resulting magnitude responses when an FIR kernel is formed using $M$ as the magnitude spectrum and one of the $P_i$ as the corresponding phase spectrum. $n$ was chosen to be $16$ here. The three magnitude responses all go through the fixed set of points defined by $M$, but vary considerably everywhere else, due to their different phase spectra.

There are many audio applications where we do not really care about the phase response of a filter, but we absolutely do care about the magnitude response. Given a magnitude spectrum $M$, we can hence regard all the $p_i$ as degrees of freedom to shape our magnitude response. More specifically, this should allow us to define $n$ additional magnitude points that the magnitude response of $K$ is required to go through. These additional requirements should be able to help to substantially reduce any undesirable magnitude response behavior of $K$ without the need to increase its length $n$. Calculating $P$ to satisfy such additional requirements seems nontrivial, though.

Are my assumptions correct? Am I missing something? Is there any practical (polynomial complexity) way to do this?

Edit:

This seems to come down to the question of wether a complex polynomial can feasibly be interpolated by a number of absolute value requirements. The coefficients of the FIR kernel $K$ define a polynomial in the complex plane and $K$'s frequency response lies on the unit circle. Any magnitude requirements are positions on the unit circle with a value which the absolute value of that polynomial is supposed to go trough. Consequently, I asked a more specific question on Mathematics Stack Exchange. If I gain any further insight from there, I will update or answer my own question here.

There are a few things that can make this easier

If you don't care about the phase response you can typically get the best result with calculating a minimum phase for your amplitude spectrum and using that as the complex target
If you want to control the in-between and the behavior at the band edges you need to specify targets points there as well.
If you have more target points then FIR coefficients, you have an over determined system but you can still solve it with a least square error approach. this will not hit all your target points exactly but it's often a good compromise between filter fidelity and computational complexity.
A least square error approach has the added advantage that you can apply weights to steer the solution towards a trade off that fits your requirements. For example, you typically need to control the out of  band behavior: you don't care about the exact behavior other than it shouldn't have too much gain. To achieve this, you can set the target simply to zero and play around with the weight until the gain in that band feels acceptable


Finding the coefficients of the block diagram I was looking at the question from a textbook of mine and it doesn't make too much sense to me:

Design a ﬁrst-order ﬁlter, linearly extrapolating the current and
previous measurement. This is equivalent to the assumption that the
upcoming measurement will be on a straight line connecting the current
and previous measurement. The output of the ﬁlter shall anticipate
this by giving this output now.

The coefficients (2 and -1) are given in the answer, but I still don't understand how they were found.
I know that the general equation of a digital filter is given as:
\begin{equation}
    y[k] = \sum_{i=0}^{N} b_i u[k-i] - \sum_{j=1}^{N} a_i y[k-j]
\end{equation}
How does one go from that, using the description of the problem, to:
\begin{equation}
    y[k] = 2 x[k] - x[k-1]
\end{equation}
This is the answer given in the book:

If you assume that the samples lie on a line then the current sample is related to the previous sample by
$$x[n]=x[n-1]+c\tag{1}$$
with some unknown constant $c$. If there are two past samples given, then $c$ is simply determined by the difference between these two samples:
$$c= x[n-1]-x[n-2]\tag{2}$$
Combining $(1)$ and $(2)$ gives
$$x[n]=2x[n-1]-x[n-2]\tag{3}$$


Generate buzzer sound I'm trying to generate a buzzer sound like this one in matlab.

When you look at the waveform, it's got sinusoidal notches cut into an otherwise normal looking sinewave.

I generate a sinewave tone,

fs = 44100 ;
t = 0 : (1/fs) : 2 ;
f = 200 ;
y = 0.5*(sin( f*2*pi*t )) ;
sound( y, fs ) ;

I don't know how to cut the notches! How would you do it, so it sounds "buzzy"?

The spectrum:

Based on the spectrum, it looks to me like this was generated by FM synthesis. The mathematics of starting with a sound and working backwards to generate it using FM synthesis can be time consuming. The best reference I know of is Computer Music. I've seen pdf editions of this book online that I thought were legal, but I'm only seeing illegal looking things now, so I'm not linking.

Some simpler approximations that might work, but almost certainly won't sound the same:

white noise with a bandpass filter. (maybe plus a little clipped sine wave)
a clipped sine wave or square wave through a low pass filter. (subtractive synthesis)


Implementation of basic 2-pole low-pass filter I know this question has been asked to hell and back from many people, and that it will probably get downvoted until it's deleted, but please bear with me.

I am designing a sound engine for one of my projects, which needs to support low-pass filtering (no other filtering required, so no need for extremely robust equalization or high-pass filters).

Unfortunately, I'm not exactly well-versed in DSP math (as I'm more a programmer than anything else), so this has been getting tricky.

What I've found by searching on the internet is a "1-RC and C filter" which is implemented as:

// Coefficient computation
// Cutoff and reso are [0,128) integers
c = 0.5^[(128-Cutoff) / 16.0]
r = 0.5^[( 24+Reso)   / 16.0]

v0 = (1.0 - r*c)*v0 - c*v1 + c*input
v1 = (1.0 - r*c)*v1 + c*v0

On the same website, I also found that c can be written frequency-wise with the formula

c = 2.0 * sin(Freq * pi/SampRate)

Since I'm only doing variable low-pass filtering with no resonance, this leads to

c = 2.0 * sin(Freq * pi/SampRate)
r = pow(0.5, 24 / 16.0) = 1 / sqrt(8) = ~0.3536

Since the sound engine is mostly MIDI-related, I've used the frequency cutoffs I've found elsewhere on the internet, which led to a relatively simple formula of

CutoffFreq = 250.0 * 32^x

Where x is bound to [0.0,1.0], which leads to a frequency range of 250Hz~8kHz.

However, this is where my issues start.

Let's assume that we're using a maximum cutoff frequency (that is, there should be no filtering on the data) and that I'm mixing at 44.1kHz.

c   = 2.0 * sin(8000.0 * pi/44100.0) = ~1.0791
r   = 8^(-0.5)                       = ~0.3536
      1.0 - r*c                      = ~0.6184

Now let's assume the input signal is a basic triangle:

Input[] = {0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 0.8, 0.6, 0.4, 0.2}

Applying the code above leads to

// Beforehand...
v0 = 0.0
v1 = 0.0

// 1st Sample [works okay]
v0 = 0.6184*0.0 - 1.0791*0.0 + 1.0791*0.0 = 0.0
v1 = 0.6184*0.0 + 1.0791*0.0              = 0.0

// 2nd Sample [incorrect]
v0 = 0.6184*0.0 - 1.0791*0.0 + 1.0791*0.2 = ~0.2158
v1 = 0.6184*0.0 + 1.0791*0.2158           = ~0.2329 [residual 0.0329]

// 3rd Sample [more than doubly incorrect]
v0 = 0.6184*0.2158 - 1.0791*0.2329 + 1.0791*0.4 = ~0.3138
v1 = 0.6184*0.2329 + 1.0791*0.3138              = ~0.4826 [residual 0.0826]

As you can see, this doesn't exactly work correctly.
In fact, judging from the sound of it, it's as though I had introduced a slight high-pass filter as well as increased resonance slightly.

As such, I was wondering if anyone could offer an explanation a more programming-related person could understand for computing low-pass filter coefficients and equations, as well as shed some light on how to properly implement this.

Thanks in advance.

First of all, your filter is a non-typical lowpass filter. The transfer function of your filter is

$H(z) = \frac{\beta^2}{1+(\beta^2-2\alpha)z^{-1}+\alpha^2z^{-2}}$

where $\alpha = 1-rc$, $\beta=c$.

In the example you use these values

$c = 2\sin(8000\pi/44100)$, $r = 8^{-0.5}$.

They give the following frequency response.

In the example where you step-by-step process the first samples of your triangular waveform it seems that you expect that input to pass through undistorted. Is that correct? if so, why do you expect this?

As an alternative you can use a 2nd order butterworth low-pass filter. Here http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt you find the expressions for the coefficients as well as the difference equation you have to implement.


Magnitude value of filter frequency response Can someone help me with this problem:

Suppose we filter signals with the below pre-emphasis filter:

$$y[n]=x[n]-0.8x[n-1]$$

I have to calculate the impulse response of the filter (ok easy) and then I have to find the magnitude of Frequency Response (dB) at $F = 2000 Hz$ with $Fs = 16000 Hz$

Can someone tell me the formula? I am confused :(

Thanks

Your filter transfer function in the Z-transform domain is given by:

$$H(z)=1-0.8z^{-1} $$

Substituting $z=e^{j\omega} $:

$$H(j\omega)=1-0.8e^{-j\omega} $$

Knowing you the sampling frequency and frequency of interest you can calculate the angular frequency: $\omega_0= 2\pi\frac{2000}{16000}=\dfrac{\pi}{4}$

Substituting you get:

$$H(j\omega_0)=1-0.8e^{-j\dfrac{\pi}{4}}$$

Taking the magnitude:

$$|H(j\omega_0)| = \left|1-0.8\dfrac{1-i}{\sqrt{2}}\right|\approx 0.7132 $$

Converting to decibel scale:

$$|H_{dB}(j\omega_0)|=20\log_{10}0.7132=-2.93\mathrm{dB} $$


How to pick a sampling frequency for Matlab fft I am recording vibration data using an adxl335 accelerometer.I have ~50ms delay before i save a reading from this sensor so most probably each second i have ~20 readings. I verified that by checking the difference between two consecutive readings since i am recording the time of reading, however since my pc is not real time i guess matlab is varying but i found that the average time between two consecutive readings is 0.0517 so my sampling should be ~19 hz ? is this the value that should be chosen ?

The sensor is analogue and being input to Arduino for analogue to digital conversion. Arduino has a function Micros() which Returns the number of microseconds since the Arduino board began running the current program, I suggest using this to accurately save the time of analogue read on Arduino. Calculate the difference as shown in the code below and send it through with the sensor reading to Matlab. Then just average the values of the time difference calculated by t2-t1 which will be almost constant, and you will have as accurate as possible the sampling time/Sampling frequency used and so you can very accurately do an FFT analysis using this information. Hope this will be helpful.

//Define Analogue Ports for Sensor input
const int xpin1 = A2; // x-axis of the accelerometer 1

//Define variable t1 for timing
unsigned long t1 = micros();

void loop()
{
  int x1 = analogRead(xpin1);  //read x from 1

  Serial.print(x1);  //print x1 value on serial monitor
  Serial.print("\t");

  unsigned long t2 = micros();
  Serial.print(t2-t1);
  Serial.print("\t");
  t1 = t2;
}


Vivado Design Suite Reed-Solomon Decoder in the description of Reed-Solomon Decoder v9.0 p 20

ERR_CNT Element
The ERR_CNT element gives the number of errors, erasures, and
punctures that were corrected. The width of the element depends on the
input parameters n and k. The width is equal to the number of binary
bits required to represent (n-k). If n-k = 16, for example, the
ERR_CNT element is five bits wide

Could someone explain me the example. 'If n-k = 16, for example, the ERR_CNT element is five bits wide'.
I thought it will be 8=(n-k)/2
EDIT 1: In the example n-k =16, then the number of errors ( error correction capability) is (n-k)/2 or I can compute it as 2*e+E<=(n-k), where e - errors, E- erasures

$\left\lfloor\frac{n-k}2\right\rfloor$ is the maximum number of bit errors that we can guarantee a code can correct for a binary channel.
But that doesn't mean a code can't sometimes correct more! (In general: not all code words have a distance that's equal to the code's minimum distance. It is the case for minimum distance separable codes, like RS codes are, though.)
Then:
RS codes are also used in applications where erasures, not bit flips, happen. You'll find the decoder documentation contains a section on erasure decoding.
You'll also find that a RS code can correct up to $n-k$ erasures – and my guess is that they made the count of that and the count of corrected errors the same width for design consistency reasons, and maybe to make different decoder IP interchangeable.
Note that giving a unsigned number more bits than it needs has no functional or FPGA utilization downside in IP cores: if whatever takes the output uses less bits than it has, the synthesizer will optimize away all calculations resulting in the bits that weren't used.


Optimizing FFT code for 12 bit processing on microcontroller I want to port the following FFT code on a microcontoller having a very little floating point capability. I want to perform FFT of the output from 12 bit ADC. Since my values are all 12 bit numbers it is wasteful of storage and processing power (I do not need very accurate FFT output) if I do every thing in double as the following code does.

How do I modify this C code so that it works only for 12 bit numbers and do not do any floating point arithmatic. This way I can save lots of storage space and  do very litlle arithmatics.

#define PI  M_PI    /* pi to machine precision, defined in math.h */
#define TWOPI   (2.0*PI)

void four1(double data[], int nn, int isign)
{
    int n, mmax, m, j, istep, i;
    double wtemp, wr, wpr, wpi, wi, theta;
    double tempr, tempi;

    n = nn << 1;
    j = 1;
    for (i = 1; i < n; i += 2) {
    if (j > i) {
        tempr = data[j];     data[j] = data[i];     data[i] = tempr;
        tempr = data[j+1]; data[j+1] = data[i+1]; data[i+1] = tempr;
    }
    m = n >> 1;
    while (m >= 2 && j > m) {
        j -= m;
        m >>= 1;
    }
    j += m;
    }
    mmax = 2;
    while (n > mmax) {
    istep = 2*mmax;
    theta = TWOPI/(isign*mmax);
    wtemp = sin(0.5*theta);
    wpr = -2.0*wtemp*wtemp;
    wpi = sin(theta);
    wr = 1.0;
    wi = 0.0;
    for (m = 1; m < mmax; m += 2) {
        for (i = m; i <= n; i += istep) {
        j =i + mmax;
        tempr = wr*data[j]   - wi*data[j+1];
        tempi = wr*data[j+1] + wi*data[j];
        data[j]   = data[i]   - tempr;
        data[j+1] = data[i+1] - tempi;
        data[i] += tempr;
        data[i+1] += tempi;
        }
        wr = (wtemp = wr)*wpr - wi*wpi + wr;
        wi = wi*wpr + wtemp*wpi + wi;
    }
    mmax = istep;
    }
}

12-bit is not a native data type on your platform, so the best you can do without loss of precision is 16-bit per sample, in Q15 format.

I don't know the Stellaris line well, but from what I gathered on the TI website, these are either Cortex M3 (LM3xxx) or Cortex M4 (LM4xxx).

The Cortex M4 embeds a floating point unit; so there won't be as much benefit in terms of run-time from converting your code to fixed point as you think. The immediate thing is to switch to 32-bit floats. The next thing is to check your compiler settings and have a quick look at the disassembled generated code to check that it actually generates floating point instructions rather than using emulation. With gcc, you'll probably have to use compilation flags like -mcpu=cortex-m4 -mfloat-abi=hard -mfpu=fpv4-sp-d16 (which are the flags I use with another Cortex-M4 product). There's a good chance you'll get better performance from using a vendor-supplied function rather than your own code. Try arm_cfft_radix4_f32 in the CMSIS DSP library. This should solve your performance problems.

If you use a M3 part instead, and/or are really short on RAM for your samples, you'll have to switch to fixed point (Q15 format). Again, you don't have to write code for that, just grab the CMSIS DSP library provided by TI and use the functions arm_rfft_q15 (for real input) or arm_cfft_radix4_q15 (for complex input, if the input size is a power of 4).